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ABSTRACT 

The structure of planetary systems around their host stars depends on their initial 
formation conditions. Massive planets will likely be formed as a consequence of rapid 
migration of planetesimals and low mass cores into specific trapping sites in proto- 
planetary discs. We present analytical modeling of inhomogeneities in protoplanetary 
discs around a variety of young stars, - from Herbig Ae/Be to classical T Tauri and 
down to M stars, - and show how they give rise to planet traps. The positions of 
these traps define the initial orbital distribution of multiple protoplanets. We investi- 
gate both corotation and Lindblad torques, and show that a new trap arises from the 
(entropy-related) corotation torque. This arises at that disc radius where disc heating 
changes from viscous to stellar irradiation dominated processes. We demonstrate that 
up to three traps (heat transitions, ice lines and dead zones) can exist in a single disc, 
and that they move differently as the disc accretion rate M decreases with time. The 
interaction between the giant planets which grow in such traps may be a crucial in- 
gredient for establishing planetary systems. We also demonstrate that the position of 
planet traps strongly depends on stellar masses and disc accretion rates. This indicates 
that host stars establish preferred scales of planetary systems formed around them. 
We discuss the potential of planet traps induced by ice lines of various molecules such 
as water and CO, and estimate the maximum and minimum mass of planets which 
undergo type I migration. We finally apply our analyses to accounting for the initial 
conditions proposed in the Nice model for the origin of our Solar system. 

Key words: accretion, accretion discs - turbulence - planets and satellites: formation 
- planet-disc interactions - protoplanetary discs - (stars:) planetary systems 



1 INTRODUCTION 

Observations of exoplanets (nearly ~ 1700 if candidates 
are included) show that the formation of multiple planets 
around their host stars is relatively commonQ This trend is 
confirmed by both the radial vel ocity and transit tec hniques 
such as the Kepler mission (e.g. [Howard et al]|201ll '). Many 
previous studies based on N-body simulation s, investigated 
the d ynamics of the planetary systems (e.g. iRasio fc FordI 
1 19961 ). It is well known that these simulations can reproduce 
the observed distribution of eccentricities of exoplanets very 
well, if they adopt a specific initial condi tion that planets are 
closely packed fe.g. iFord fc Rasidll2008l ). For our Solar sys- 
tem, the Nice model which requires a specific initial arrange- 
ment of Jupiter and Saturn explains the dynamics very well 
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(|Morbidellill2010l . references herein). The obvious question is 
whether or not such conditions are reasonable. In order to 
answer the question, one needs to consider the early stage of 
planet formation in which phy sical processes ar e controlled 
by gas dynamics in discs (e.g. iThommes et al.l 1^008). The 
presence of gas is necessary for gas giants to be formed. 
However, it also causes rapid inward planetary migration 
(|Wardlll997l '). Since this radial motion strongly depends on 
the disc properties such as the gas surface density and the 
disc temperature, it becomes a huge challenge to systemat- 
ically investigate what are the realistic initial conditions for 
the later evolution of planetary systems. 

Rece ntly, planet traps have received a lot of atten- 
tion (M asset et all l2006l : iMorbidelli et all ,2008). Barriers 
to planetary migration arise because of inhomogeneities in 
discs where the direction of planetary migration switches 
from inwards to outwards, so that migrating planets are 
halted. (Equivalently, the net torque exerted on planets is 
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zero there, iHasegawa fc PudritjIioiOal . hereafter HPIO, ref- 
erences herein). Planets may acquire most of their mass as 
they accrete material at these barriers - which we call planet 
traps in globally evolving discs. We use the term barriers for 
our local analyses (§[3]and|4ll while the term of planet traps 
are used for our global, unified analyses (§0- 

Planet traps were originally proposed bv lMasset et al.l 

l|2006t) in order to solve the well known rapid migration prob- 
lem wherein plan e ts can be lost to discs within 10^ years. 
iMatsumura et all l|2007l . hereafter MPT07) first addressed 
a (possible) link between the planet traps and the diver- 
sity of exoplanets by showing that planet traps can move 
du e to the time dependent , viscous evolution of discs (also 
see IMatsumura et al.ll2009l ). Although they focused on dead 
zones in discs (which are the high density, inner regions, so 
that turbulence there induced by magnetorotational insta- 
bility (MRI) is quenched, ICammic 1996), their basic idea is 
valid for any type of planet trap. More recently, it has been 
pointe d out that a sin g le disc can have a couple of planet 
traps l|lda fc LinI l2008l : iLvra et all l2O10h . Population syn- 
thesis models confirm that the planet traps and their move- 
ments can be imp ortant for the diversit y of observed exo- 
planetary systems (|Mordasini et al.l2011^ . As already noted, 
planet traps have the possibilit y of strongly enhan cing the 
growth rate of planetary cores jSandor et al.|[201ll ) and the 
formation of giant planets. Thus, planet traps have consider- 
able potential for understanding the formation of planetary 
systems. 

In a series of papers, we systematically investigate how 
inhomogeneities in discs can trap migrating planets in over- 
dense regions where they undergo most of their growth, and 
how trapped planets establish their planetary systems in 
viscously evolving discs. In this paper, we undertake a com- 
prehensive study of the various mechanisms that produce 
planet traps in discs around a variety of young stars, - from 
high to low mass. One of our major findings is that a single 
disc can have up to three planet trap regions. We show that 
the position of any trap depends upon the disc's accretion 
rate, which decreases with time as the disc is used up and 
star formation is terminated. The decreasing accretion rate 
forces the planet traps to move inwards at varying rates. 
This ultimately sets up the condition for the mutual inter- 
action of planets in the traps, which can provide the realistic 
initial conditions for the evolution of planetary systems. 

The plan of this paper is as follows. We describe physi- 
cal processes that create inhomogeneities in discs and sum- 
marise our analytical approach for them in § [2] The general 
reader may then proceed to §[S]for discussion of the general 
results (also see Fig. [7]). Armed with physical understand- 
ing of the inhomogeneities, we investigate corotation and 
Lindblad torques which are the driving force of planetary 
migration in §[3]and|4l respectively. For the former case, we 
show that a new barrier arises due to the change of the main 
heat source for gas disc from viscosity to stellar irradiation. 
For both torques , ice lines play s ome role in trapp ing planets 
l|lda fc Linll2008l . hereafter IL08: lLvra et al.ll2010l ). In §[3 we 
discuss other possible barriers and estimate the maximum 
and minimum mass with which planets undergo type I mi- 
gration. Also, we investigate the possibility of the presence 
of ice lines due to various molecules. In § [SJ we integrate 
our analyses and discuss the roles of these planet traps in 
the formation of planetary systems. We apply our results to 



an explanation of the Nice model for the architecture of the 
Solar system in § [T] In § [S] we present our conclusions. We 
summarise important quantities which often appear in this 
paper in Table [T] 



2 DISC INHOMOGENEITIES 

We describe physical processes governing the structure of 
protoplanetary discs and discuss how disc inhomogeneities 
arises from these processes. We discuss the basic features 
of our analytical modeling of the resultant disc structures - 
which will be presented in § |3] and [l] Table [2] summarises 
the disc inhomogeneities, related nomenclature that often 
appears in the literature, the dominant torque that trans- 
ports angular momentum in that region of the disc, and the 
section of the paper that treats the analysis. 

2.1 Physical processes 

Heating of protoplanetary discs is one of the most important 
proc esses in order to unders tand their geometrical structure 
(e.g. iDuUemond et al.ll2007h . Since discs are accreted onto 
the central stars, release of gravitational energy through disc 
accretion becomes one of the main heat sources. This en- 
ergy can be dissipated by viscous stresses, leading to vis- 
cous heating. Once discs are heated up, then the absorption 
efficiency of discs which is regulated by their optical depth 
establishes the thermal structure of discs. In protoplanetary 
discs, dust gives the main contribution to opacity for pho- 
tons with low to intermediate energy while gas is the main 
absorber of high energy photons. In the inner region of discs, 
viscous heating is very efficient and leads to high disc tem- 
peratures (> 1000 K for classical T Tauri stars (CTTSs), 
iD'Alessio et al.lll998l ). As a result, both metal dust grain s 
and molecules are destroyed there (e.g. iBeU fc LinI 1 19941 ). 
Thus, this high disc temperature and resultant destruction 
of opacity sources produces opacity transitions. It is inter- 
esting that opacity transitions are also produced in a com- 
pletely different situation, that is, of low disc temperatures 
wherein opacity is enhanced by "freeze-out" processes. Ice 
lines are one of the most famous opacity transitions created 
by this process. Ice lines are defined such that the number 
density of ice-coated dust grains suddenly increases, which 
is a consequence of low disc temperatures there. In § 13.2.41 
we treat ice lines as an example of opacity transitions. 

It is well known that protoplanetar y discs are also 
heated up by stellar irradia tion (Chiang fc Goldreichlll997l : 
iHasegawa fc Pudritj[2010bh . Since viscous heating becomes 
less efficient for larger disc radii, stellar irradiation plays a 
dominant role in regulating the thermal structure of discs 
there. Thus, a transition exists for discs wherein the dom- 
inant heating mechanism transits from viscous heating to 
stellar irradiation. We call this a heat transition, (see § 

[323}. 

We have focused so far on heating processes which in- 
volve with photons that have relatively low energy. High 
energy photons such as X-rays from the central stars and 
cosmic rays, are important for understanding the ionization 
of protoplanetary discs. Since the MRI is the most favoured 
process to excite turbulence in discs , it is crucial t o evalu- 
ate the ionization structure of discs (jGammiell 19961 ). In the 
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Table 1. Important quantities 



Symbols 


Meaning 


A/* 


Stellar ma^s 


R* 


Stellar radius 


T* 


Stellar effective temperature 


M 


Accretion rate (see ec^uation (j3.6[l ) 


G 


Gravitational constant 


n 


Angular frequency (see equation ||4.3[l) 


K ep 


Keplerian frequency (= ■\JG]\''I^It'~'^ 


IVlp 


Planetary mass 


p 


Planetary orbital radius 


11 




n 
r 


Gas volume density 


Si 


Gas surface density 








Power-law index of S if S oc 


T 


Disc temperature (cx! r'') 


t 


Power- law index of T 


T 






Disc temperature of the mid-plane 




Condensation temperature for species k at the ice line 


Ts 


Disc temperature of the surface 


Cs 




H 




H 


Disc photosphere height 


h 


Disc aspect ratio (^ Hjv^ 


'^trans 


Disc radius of disc inhomogeneities 


Fr. 


Gas density modification at 'i^trans 




Transition width at vtvans cH'^ 




Disc radius of the heat transition 




Disc radius of ice lines 


'^edge 


Disc radius of the outer edge of dead zones 


P 


(^nc; TM*fiGGi Iff ( — r>f'^\ 

VjrcLD J^l CDO Li. J. " ^ fjl^ Q J 




Epicyclic frequency (see equation (]4.4[1 J 




V V CL V dl Lllil L/Ci. loCC CLJ LiCLUlWll Jrr i ^ |^ J 




T,1Tl/Hr^lf^*H T*fCJO"r\QTl'i" T^/^CI'i'l/^Tl i 1 1 

J-JlimUltlU. 1 COUlldllL JptJOl LrltJli 1 p ) 


*r 


Forcins function fsee eo nations (14. 5Il and (14. 12 II ") 




V dl Lit; Ul CL L^LldliLlL^ v\. CLL / p 




kJ Lll. IfXL^^ LJ-L^llolLV Wl dL. LI Vc Iti^lLJllO lOcC CLJLtdLlLJll l|rc*^rrUf 


^ A 


Power-law index of (generally ^ 0) 


ro 


V_/lltlyl dl^LCl iO LlLj LllBLj 1 dLli LIB 


f. 

J ice 


A/TrtH 1 fir"a 1 1 rtvi rnf+OT* ot t! ^ Hup to iff Iit^pg 1 — tIt' i r n r iiill 

lV±LJLi.IllLjd LlLJll IdL- LLJl LlLiC LLJ IL^C llllCo I / \l i[ ^ J dl. ^ J d2, } i 


J dl 


Parameter representing density jumps due to ice lines 


J dl 


Parameter representing dust traps due to ice lines 


9d 


Parameter representing density bumps due to ice lines 


a 


Mean strength of turbulence fsee equation (14.21Il) 


OLA 


Strength of turbulence in the active zone (= 10~^) 


OLD 


Strength of turbulence in the dead zone (= 10~^) 


V 


Kinematic viscosity (= aCgH) 


7 


Adiabatic index (=1.4) 




Heat flux in the vertical direction 




Viscous dissipation rate per unit mass 


O-SB 


the Stefan-Boltzmann constant 


Ub 


the Boltzmann constant 


(^GA 


Grazing angle fsee equation (13.20ll) 


R 


Opacity 


T 


Optical depth (= kS) 



1 The simplest assumption that S oc r** is only adopted in I3.2l and l4.2l 



inner region of discs, the column density is very high, so 
that high energy photons cannot penetrate the entire re- 
gion. In this "layered" structure, the mid-plane is effectively 
screened from being ionized while the surface region is effi- 
ciently ionized. Thus, the inner region is less turbulent, and 



characterised by a dead zone. On the other hand, the outer 
region has lower column density and hence is fully ionized 
for the entire region including the mid-plane. As a result, 
the outer region is fully turbulent. Thus, transitions in the 
amplitude of turbulence exist in discs which we denote as 
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Table 2. Summary of disc inhomogeneities 



Disc inhomogcneity 


Nomenclature 


Torque 


Section 


Opacity transitions 


e.g. Ice lines ^ 


Corotation 


13.2.41 


Heat transitions 


N/A 


Corotation 


|3.2.3| 


Turbulence transitions 


Dead zones 


Lindblad 


1431 


Opacity and turbulent transitions 


Ice lines^ 


Lindbald 


m 



^ Disc turbulence is assumed to originate from non-MRI based processes 
^ Disc turbulence is assumed to originate from MRI based processes 



turbulent transitions (see § 14. 3p . Dead zones are the most 
famous product of the turbulent transitions. 

Finally, we introduce interesting transitions which arise 
from the combination of opacity and turbulent transitions. If 
the MRI is the mechanism that drives disc turbulence, then 
ice lines can also be regarded as this type of transition (IL08, 
see 14.4. II for the complete discussion). At the ice lines, the 
sudden increment of ice-coated dust grains is likely to result 
in the sudden decrement of free electrons there, since such 
sticky grains can efficiently absorb them (fSano et al.ll200ol ). 
This process, initiated by the opacity transitions, reduces 
the ionization level there, and therefore removes the coupling 
between the magnetic field and the gas - which kills the MRI 
instability. Hence ice lines create turbulent transitions. In 
summary, ice lines are regarded as an opacity transition if 
disc turbulence is excited by non MRI processes while they 
act as an opacity and turbulent transition if turbulence is 
excited by the MRI. In § 14.41 we investigate ice lines again, 
but as an example of the opacity and turbulent transitions. 



3 BARRIERS ARISING FROM COROTATION 
TORQUES 

Corotation torques play an important role in planetary mi- 
gration in regions with high viscosity (10~^ ^ ^ 10~^), be- 
cause the corotation torque ther e is (partially) unsaturated 
(i.e. effective, see Appendix B in Hase gawa fc Pudrit j|201ll . 
references herein). The corotation torque is known to act as a 
barrier for t wo situations. The first situation is when S oc r" 
with s > 1 l|Masset et al.|[200^ : iPaardekooper fc Papaloizoul 
This case can be generally established only at the in- 
ner edge of discs. In typical discs, the inner edge is located 
around a few hundredths of an, which corresponds to the 
semi-major axis of Hot Jupiters. Thus, the positive gradi- 
ents of surface density unlikely play a dominant role in ex- 
plaining the diversity of the observed planetary systems. The 
second sit uation arises for adiabatic discs that have inhomo- 
geneities (|Lvra et al.ll201Q ). This is more promising since the 
main sites of planetary formation in protoplanetary discs are 
generally optically thick and reasonably approximated to be 
adiabatic. In this section, we investigate only the second case 
for the above reason. 



3.1 Conditions for outward migration 



2.2 Our analytical approach 

We present analytical modeling of the disc structures af- 
fected by the disc inhomogeneities in § [3] and |4l In order to 
make analytical treatments possible, we adopt expressions 
that are simple enough, but well capture the physics aris- 
ing from the inhomogeneities. As an example, we adopt a 
tanh function for representing an inhomogcneity in the sur- 
face density (see equation (I3.1[l ). This profile is more general 
than a power-law and is likely to be applicable for the cases 
of opacity, heat, and turbulent transi tions. For the opac- 
ity an d heat transitions, the results of iMenou fc GoodmanI 
l|2004l . hereafter MG04, see their figl) validate the usage of 
this function as do the results of MPT07 for the turbulent 
transitions. Disc structures are approximated as power-laws 
for regions far away from the inhomogeneities (see equation 
((ST}). Thus, we adopt more general profiles for the disc 
structures. In addition, we make use of the simplest expres- 
sions in order to reduce mathematical complexity, and hence 
different functions are adopted for different transitions. 

Bearing these in mind, we will demonstrate that the 
disc inhomogeneities are the most plausible sites to produce 
barriers to type I migration by undertaking a comprehensive 
study of both corotation and Lindblad torques in §[3]and|4l 
respectively. 



We derive conditions for the disc structures that are required 
for planets to migrate outwards. These conditions are crucial 
for the subsequent subsection where we discuss how disc 
inhomogeneities work as a barrier. 

3.1.1 Disc models 

We discuss our disc models (also see Table [T|. For the sur- 
face density affected by disc inhomogeneities, we adopt the 
following equation; 



1 + 



Fo 



1 



(l-tanh(!l^)) 



(3.1) 



where Si„i oc r" is the initial density profile, Fg charac- 
terise the density distortion which is a consequence of disc 
inhomogeneities, rtrans is a orbital radius of the disc inho- 
mogeneities, and uj = cH is the width of the transition (also 
see Table [T|. This analytical modeling is motivated by the 
results of MG04 who first found the significant effects of 
opacity transitions on planetary migration by solving the 
detailed, ID disc structure equations. Based on their results 
(see their fig. 1) , any temperature distortion created by opac- 
ity or heat transitions, is likely to result in a surface density 
structure which is well expressed by equation (|3.ip . For the 
initial profile (s), we examine two cases both of which are 
well discussed in the literature: s = —3/2 and s — —1. The 
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former power-law index is known as minimum mass solar 
nebula (MMSN) while the latter one is a steady state solu- 
tion to disc accretion. For the disc temperature, we adopt 
simple, power-law structures (T oc r') and derive the critical 
value of t that results in outward migration below. 



3.1.2 Torque formula 

We ad opt the torque formula derived bv lPaardekooper et al.l 
l|20ld ) in which discs are assumed to be 2D. In the for- 
mula, the total torque is comprised of the linear Lindblad 
torque and non-linear corotation torques, known as horse- 
shoe drags. If discs are (locally) isothermal, both Lind- 
blad and corotation torques dictate that standard rapid in- 
ward migration will occur. More specifically, the vortensity- 
related horseshoe drag that is only an active (corotation) 
torque in this case , cannot exceed the Lindblad t orque 
for power-law discs (|Paardekooper &: Papaloizoull2009l '). For 
discs with more general profiles, the vortensity-related horse- 
shoe drag can exceed the Lindblad torque, but still re- 
sults in inward migration. Therefore, we assume discs to 
be adiabatic, which is reasonable fo r the inner part of 
discs (r < 80 an) around CTTSs (|Chiang fc GoldreichI 
19971). Adopting the torq ue formula in adiabatic discs (see 
Paardekooper et al.ll201Cl ). the direction of migration may 



be given as 



sgn 



(2.5 - l.lt + Q.ls) 



1.10, {l+s^+I^(t-{^-l)s) 

^ (3.2) " 

where s = dlnE/dlnr is an "effective" power-law index 
for E which can have a more general profile, 7 is the adia- 
batic index and we took the softening length h = Q.Ah for 
simplicity (also see Table [TJ. The terms in the first brack- 
ets arise from the Lindblad torque, the terms in the second 
from the vortensity-related horseshoe drag, and the terms 
in the third from the en t ropy-r elated horseshoe drag. As 
iPaardekooper fc Mellemal (|2006l ) first discovered and sub- 
sequent works interpreted it, the entropy-related horseshoe 
drag which occurs only in adiabatic discs, is scaled by ra- 
dial, entropy gradients and c an drive planets into outw ard 
migration (see Appendix B in lHasegawa fc Pudrit j[201ll for 

a summary, references herein). 

The original formula derived by IPaardekooper et al.l 

is valid exclusively in power-law discs. In order to 
take into account discs with more general profiles such as 
equation (|3.1|l . we add a vortensity correction factor 01, in 
the second term which is defined as 



dln(S/B) 



where 



s-f3/2 dlnr 



5 = 1^(^2^) 
2r dr^ ' 



(3.3) 



(3.4) 



is the vorticity of gas, also known as one of the 
Oort constants. This is because gradients of vortensities 
{d\n{Yi/ B) / dlnr) are very sensitive to disc structure0 We 



^ Lindblad torques also need a similar correction factor for discs 
with general profiles. Indeed, we derive an analytical relation 
for Lindblad torques using equation 13.11 1 in § 14.31 (see equation 



note that = 1 for pure power-law behaviours. In fact, gra- 
dients of vortensities are the core of the vortensity-related 
corotation torque and regulate the transf er of angular mo- 
ment um there (see Appendix B in iHasegawa fc Pudritd 
I2OIII ). Therefore, inclusion of the factor 0„ is likely to be cru- 
cial for properly evaluating the importance of the vortensity- 
related corotation torque relative to the others. 

The condition for outward migration is, therefore. 



t < 



[1.10,(3/2 + s)) - 2.5 + 7.8s]7 - 7.9s 
-1.774-7.9 



(3.5) 



It is fruitful to first consider pure power-law discs where 
01, = 1. Setting 7 = 1.4, the required temperature profile is 
t < -1.4 for MMSN discs (s = -3/2) while for s = -1, the 
temperature profile t < — 1.1 is needed. We stress that such 
steep temperature profiles can be only achieved by viscous 
heating in optically thick discs, which is discussed more in 
the next subsection. 

Let us now investigate how the required t deviates from 
the power-law predictions due to the factor 0,. In general, 
we find that the structure of 0„ is very complicated. There- 
fore, we present the detail discussion of 0„ in Appendix [K\ 
and briefiy summarise the two most important effects on 
t here. The first is that the effects of 0i, are very local 
and are likely to be confined within the transition region 
A r ~ a; = cH. This is clearly shown in Fig. [T] (see the ver- 
tical dotted line on the bottom panel for representing the 
transition region). In this figure, we set c = 1 that is the 
most likely value for the opacity and heat transitions. For 
the bottom panel, we plot the critical value of t derived from 
equation (|3.5|) as a function of rtrans/ro, where the charac- 
teristic disc radius is set ro = 1 au unless otherwise stated. 
The top panel shows the behaviour of s. For both panels, the 
solid line denotes the results for the disc with the general 
profiles (see equation (|3.1[) ) while the dashed line is for the 
pure power-law discs. We show the results only for the case 
that T,i„t oc r~^, because the results of Ei„t oc r~^^^ are 
qualitatively similar. The second important point is that, 
in the local region where the effects of 0„ become crucial, 
a required t reaches values that are unattainable by any 
physical process in protoplanetary discs. This means that 
outward migration cannot happen due to large vortensity- 
related corotation torques that results in inward migration. 
In summary, the inclusion of 0„ reduces the possible region 
wherein planets can migrate outwards, but this reduction is 
well confined in a local region (~ an order of lj) centered at 

r — rtrans. 

Finally, we note that this torque formula does not take 
into account any effect of saturation. The problem of satura- 
tion, which is central to the problem of corotation torque, is 
very complicated and st rongly dependent on th e fiow pattern 
at the hors eshoe orb it ('Masset fc Casolil 2010l: also see Ap- 
pendi x B in Hasegaw a & Pudritz 2011lT |Paardekooper et aP 
(j201ll ): [Masset fc Casolil ((2010) attempted to derive analyt- 
ical formulae in which the effects of saturation are included. 



I|4.20|l ). As shown in AppendixlA] however, they switch the direc- 
tion of migration only if Fg > 1.5 for c = 1. In general, the density 
distortion created by opacity and heat transitions is likely to be 
less than that. In addition, it is more consistent to use the above 
terms than equation l|4.20|l in this torque formulation. Therefore, 
wc adopt the original form for the Lindblad torques. 
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Figure 1. The effects of the vortensity correction factor 0i, on s 
and t on the top and bottom panel, respectively. For both panels, 
the solid lines denote the results for discs with general profiles 
(see equation lj3.1|l ) while the dashed lines are for pure power-law 
discs. Only the results of Si„t <x are shown because those of 
the MMSN case are qualitatively similar. We set the transition 
width parameter c = 1. For both panels, the effects of are well 
confined in the transition region oj) (see the vertical dotted 
lines on the bottom panel). In addition, inclusion of ij)^ reduces 
the possibility of outward migration there. 



However, both formulae depend on thermal diffusivity in 
discs that is totally unknown in protoplanetary discs. Hence, 
we adopted the unsaturated torque formula. This implies 
that the above required temperature profiles may be the 
minimum value. Steeper profiles may be needed if (partial) 
saturation effects are taken into account. 



3.2 Heat and opacity transition barriers 

The entropy-related horseshoe drag can result in outward 
migration only in discs with steep temperature profiles, as 
discussed above. Only viscous heating can establish such 
profiles. On the other hand, stellar irradiation can heat up 
protoplanetary discs as well. By deriving temperature pro- 
files for both heating processes, we show that the heating 
transition from viscosity to stellar irradiation activates a 
new barrier. In addition, we examine the effects of opacity 
transitions on the temperature structure. We exclusively fo- 
cus on ice lines as an example of the opacity transitions. For 
both transitions, we identify the positions of barriers. 



Table 3. Typical quantities for stars with various masses 



Her big Ae/Be stars CTTSs M stars 



M (Mq/ year) 


10-'' 


10-8 


10-10 


M. (M@) 


2.5 


0.5 


0.1 


R* (Re) 


2 


2.5 


0.4 


n (K) 


10000 


4000 


2850 


(g cm =') 


10^ 


10^ 


102 



1 The value of Eq is 10 times larger than the standard values. 



3.2.1 Disc models 

We adopt simple, power-law profiles for both the surface 
density and disc temperature (E oc r" and T oc r*). This is 
supported by the argument done in the above subsection. 
Since the effects of 4>v{^ 1) occurs only in a local region 
centered at the position of a transition, none of our findings 
discussed below is affected. In addition, the resultant devi- 
ation for our estimate of the position of barriers is only < 
10 per cents. Thus, it is reasonable to use power-law discs 
here. 

We take the value of surface density at ro (labeled by 
Eo) as 10 times larger than the standard values of discs for 
all cases. Many previous studies showed that gas giants can 
be formed only if the disc mass i ncreases by this am ount 
relative to the MMSN models (e.g. lAlibert et al.|[2005l ). Ta- 
ble [3] tabulates the values of Eo around stars with various 
masses. 

In addition, we assume stationary accretion discs that 
have accretion rates which can be modeled as 

M = Stti^E = SnaCsHT. (3.6) 

using the famous a- prescription (jShakura fc SunvaevlllOTSi ). 



3.2.2 Viscous heating: Origin of outward migration 

We first discuss viscous heating. Assuming local thermal 
equilibrium in geometrically thin discs, the energy equation 
can be written as 



p dz 



(3.7) 



where is the heat flux in the vertical direction and D^is 
is the viscous ene rgy dissipation rate per unit mass (e.g. 
iRuden fc Linlll986l . also see Table [l}. In Keplerian discs, the 
viscous dissipation rate becomes 



dQ. 
dr 



Integrating equation (|3.7|l gives 

9 



(3.8) 



(3.9) 



If discs are assumed to be optically thin and isothermal, 
the flux Fz = 2(TsBTg with the Stefan-Boltzmann constant 
crsB. Thus, the resultant effective temperature profile be- 
comes 

oc r"^/* (3.10) 
for discs which are assumed to reach a steady state. In this 
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case which is valid in the late stage of disc evolution, ac- 
cretion rates (Af) is constant over t he entire disc. This fa- 
mous 3/4 law was first derived by iLvnden-Bell &: Pringld 
and is exactly identical to the temperature pro- 
file for flat discs which are heated by stella r irradiation 
(|Adams et al.lll987l : IChiang fc Goldreich|[l997l ) . In addition, 
it is well known that the resultant spectral energy distribu- 
tions (SEDs) do not reproduce the observed ones - flared 
discs are required. 

If the accretion rate (M) changes with disc radius, then 
the self-consistent temperature proflle is 



oc r 



s/3-1/2 



(3.11) 



where E cx r° is the surface density of discs. For the MMSN 
models (s — —3/2), Te oc while for discs with s = —1, 
Te oc r~^^^ . It is obvious that the temperature profiles de- 
rived from the isothermal assumption are not steep enough 
for the corotation torque to provide a barrier. 

The isothermal assumption can break down, especially 
in the main site of planetary formation (lau < r < 20au). 
The region of discs is reasonably considered to be optically 
thick and the energy is mainly transported by radiation 
l|D'Alessio et al.lll998l ). In this case, the flux is generally de- 
scribed by 



(3.12) 



3kp dz ' 

where R is the opacity (e.g. lRuden fc Linll986l . also see Table 
[!}. Assuming the flux Fz is constant over z, integration of 
equation (|3.12[) gives the relation between the temperatures 
of surface and midplane regions; 



3r 



(3.13) 



where r = kE is the optical depth and Ts and Tm is the tem- 
per ature of surface and midplane region, respectively (also 
see iNakamoto fc Nakagawal 1 1 9941 . Table [T]) . If only viscous 
heating is taken into account, the temperature of the surface 
region should be much smaller than that of the midplane. 
As a result, the temperature structure is governed by 



-*- m. 



640-SB 



(3.14) 



When radiative t ransfer equation is treated e xplicitly, r is 
replaced by Te// (|Hubenvlll990l : iKlev fc CridalBoOS l : 



Teff = -T - 



1 

4 4r' 



If discs have a constant M (see equation (|3.6|l 
the temperature profile becomes 

rrn s/4-3/4 



(3.15) 



then 



(3.16) 



assuming R to be independent of E and T. For the MMSN 
models (s = —3/2), Tm oc r~^^^ while for discs with s — —1, 
Tm OC . Again, these temperature profiles are not steep 

enough. 

However, iKlev fc Cridal (|2008D performed numerical 
simulations by solving a more complicated energy equation, 
and showed that the resultant temperature profile goes to 
r~^'^ for discs with E oc r~^'^ in st eady state. We ca n gain a 
similar profile if we adopt R oc (|BeU fc Linll 19941 ). In this 
case, Tm oc r^^^"^/^. Therefore, the assumption that R is 



independent of E and T, likely underestimates the tempera- 
ture structure. If M is not constant, then the self-consistent 
temperature profile becomes 



Tm oc r 



2s/3-l/2 



(3.17) 



assuming R to be independent of E and T. For the MMSN 
models (s — —3/2), Tm oc r~^'^ while for discs with s = 
— 1, Tm oc r~^'^ . These profiles are exactly what is required 
(equation (|3.2|) ') in order for planets to migrate outwards. If 
K oc T^ is adopted, Tm oc r^""^/^, which is more preferred 
for outward migration. 

Thus, outward migration due to the entropy-related 
horseshoe drag is expected for viscously heated, optically 
thick discs. If protoplanetary discs were homogeneous in 
opacity and viscosity were the only physical process heating 
them, the results would suggest that large region of discs 
ought to be devoid of planets. This does not occur, however. 
We examine two kinds of inhomogeneity of discs below. One 
of them is stellar irradiation, which gives a new barrier. The 
other is the ice li ne which arises as a consequence of an 
opacity transition (iLvra et al.ll2010h . 



3.2.3 Stellar irradiation: Origin of heat transition barriers 

Stellar irradiation is well known to be the main heat source 
for regions in the di sc beyond r > 2 — 3 au (in CTTSs, 
iD'Alessio et al.lll998l ). Inside of that radius, viscous heat- 
ing dominates over stellar irradiation. Since the temperature 
slope controlled by stellar irradiation is much shallower than 
that of viscous heating, planets which migrate inwards can 
be halted at the turning point where viscous heating begins 
to take over. We call this barrier a heat transition barrier. 

In order to calculate the position of the turning point, 
we ad opt radiative disc models of IChiang fc GoldreichI 
(I1997I ). In this model, two kinds of disc temperature are 
calculated. One of them represents the surface layer which 
is directly heated by the central star, called the super-heated 
layer. The other is for the midplane layer which is heated 
by the super-heated layer. Assuming radiative equilibrium 
discs, the temperature of the super-heated layer is given as 



Ts 



so (^) 

V au/ 



-2/5 



(3.18) 



where Tso ~ 550 K for discs around the CTTSs. Since the 
super-heated layer radiates equal amounts of energy inwards 
and outwards, the temperature of the midplane layer which 
is approximated to be optically thick is written as 



V 4 



1/2 



(3.19) 



where acA is the grazing angle at which photons emitted 
from the star strike the discs (also see Table [TJ . In general, 
acA is written as 



0.4i?, d /H 

OtGA ~ h r— — 

r dr \ r 



(3.20) 



where H is the height of the visible photosphere above 
the mid-plane. Thus, the temperature of the midplane is 
strongly dependent on the geometrical shape of discs if stel- 
lar irradiation is taken into account. 

For flat discs in which the aspect ratio h is the constant. 
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Table 4. The typical values of r^jt for So 







Herbig Ae/Be stars 


CTTSs 


M stars 


r^t (au) for s 


= -3/2 


159 


16 


1.8 


rht (au) for 


s = -1 


827 


40 


2.3 



the temperature of the midplane layer becomes 

, 1/4 



(3.21) 



Again, this is the famous 3/4 law. For flared discs 
which are needed for rep roducing the observed SEDs 
l|Kenvon fc HartmannllT987l ). the second term of the right 
hand side in equation (|3.20p becomes dominant. Assuming 
that dust is well mixed with the gas in discs, 

1/2 / „ ^^ 1/2 



H 

r 



HH 



H_ 
H 



(3.22) 



where 



(3.23) 



/ig is the mean molecular weight of the gas and fc_g is the 
Boltzmann constant (also see Table [ij. Assuming H/H is 
constant, the self-consistent temperature of the midplane is 



T ( ^* 

' J-mO 



3/7 



where 



\im) 



1/7 



(3.24) 



(3.25) 



It is obvious that this temperature profile is much shallower 
than that due to viscous heating. 

Thus, the radius at which temperature profile switches 
is (equating equations (|3.14p and (|3.24l) ) 



ro 



TmO I 



3/7 



64asB^J■g 
27KoEQa7fcsno 



-3/2+3/7 



(3.26) 



where we adopt k — kq with Ko = 2 X 10"'' (jBell fc LinI 
1 19941 ). The usage of this form of R that is only valid for the 
region outside of ice lines is reasonable, which is discussed 
below. 

Fig. [2] shows the heat transition radius for discs around 
stars with various masses. For the solid lines, we adopt the 
values of Eq in Table [3] with a — a a (also see Table [T]). 
For comparison purposes, the dashed lines denote the case 
with 0.1 X Eo (which is the standard surface density in the 
literature). The transition radius generally increases with 
increasing s. Otherwise, it becomes a decreasing function 
of s (see the dashed line for the case of M stars). Table [4] 
summarises the values of rht for both s = —3/2 and s = 
— 1 with Eq. It is obvious that massive discs around Herbig 
Ae/Be stars extend Vrh to the order of several hundred au 
while lowest mass discs around M stars shrink Vht to the 
order of 1 au. For the intermediate disc mass around CTTSs, 
rht is the order of a few ten au. 



1000 
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1 00 



< 
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r^^j with Eq 

r^t with 0.1 xEq 
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M stars 



rp,( with Eq 
ri,, with 0. 



— I .D 
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.3 -1.2 - 
s 



-1 .0 



Figure 2. The heat transition radius as a function of the power- 
law index s. For aU panels, the solid line denotes the case using 
the values in Table |3] while only Sq decreases by a factor of 10 
for the dashed line, which is the value used more often in the 
literature. In general, the radius increase with increasing s, since 
the quantity in the brackets of equation I I3.26I I is larger than unity. 
If it is less than unity (see the dashed line for the case of M stars) , 
the radius decreases. 
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3.2.4 ^ce lines: Origin of ice line barriers 

We now discuss the effects of ice lin es on temperature pro- 
files which provide another barrier (jLvra et al.l l201(il '). We 
note that our following approach is applicable to any other 
opacity transitions. We discuss them more in i) l5.2l At the ice 
line, the number density of icy dust grains suddenly increases 
due to the low disc temperature. This sudden increment of 
dust provides an opacity transition. As a result, the temper- 
ature profile around the ice lines becomes shallow, resulting 
in inward migration while the profile well inside of the ice 
line becomes steep enough for planets to migrate outwards. 
Thus, planets migrating inwards are halted around the ice 
lines. This is known as the ice line barrier. This assumes 
that corotation torques are active -i.e. non-saturated, which 
is not necessarily clear at the ice lines (see § 14.4.11 and 15. 2p . 

At first, we examine which heat source controls the lo- 
cation of ice lines, ru. In order to proceed, we adopt equation 
(|3.14|) for viscous heating and equation (|3.24[) for stellar irra- 
diation. Equating Tm to the cond ensation temperature for 
H2O ice, T^;}i^o{ru) = 170 K l|jang-Condeli fc Sasselovl 
|2004| '1. equation (|3.14|) becomes 



ro 



T 



, ^ 64o-gfl/j3 
>.H20i'^«''27KoE2Q7fcsno 



(3.27) 



where k = koT with kq 
equation (|3.24|l becomes 



2 X 10^'' ijBell fc Linll 19941 1. and 



7/3 



(3.28) 



In general, the exponent index s is negative, so that ru oc 
j^2/(3/2-2|s|) ^ This shows that as the surface density (So) 
decreases, the position of the ice line barriers moves inwards. 
We note that we examine the effects of a water-ice line here 
although the above argument is applicable to ice lines of any 
material by changing the relevant condensation temperature 
for species k, Tm, *;(?■«)■ We will discuss them more in 15.31 

Fig. [3] shows the above two equations as a function of a 
for discs around stars with various masses. We set s = —1.5, 
since this choice of s minimises the importance of viscous 
heating in discs (see Fig. [2|. The black, solid and dashed 
lines denote equation (|3.27p with Eo and 0.1 x So, respec- 
tively while the dotted line is for equation (|3.28p . In all cases, 
ru defined by viscous heating is located at greater distances 
from the star than that by stellar irradiation. Furthermore, 
the transition radius defined by equation (|3.26|) is larger than 
any of these two radii (see the gray, thick lines). Thus, we 
can conclude that ru is determined by viscous heating for 
discs wi th a wide r ange of a. This agrees with numerical 
work bv lMin et al.l (|201ll ) who showed the same results by 
numerically solving the full wavelength dependent, radia- 
tive transfer equation by means of a Monte Carlo method 
in 3D discs. In their simulations, viscous heating is explic- 
itly included as well as s tellar irrad i ation. In addition, this 
supports the findings of iLvra et al.l l|2010l 'l who adopt disc 
models which are only heated by viscous heating. Further- 
more, Fig. shows that the radius of the heat transition is 
located outside of the ice line. This validates our choice of R 
in equation (|3.26|) . 

We now investigate temperature profiles created by vis- 
cous heating around ice lines by adopting appropriate R for 
each region (|Bell fc Linlll994l '). For the region inside of the 



ice lines, k — kqT^^^ with kq = 0.1. As a result, equation 
1)3. 14[) becomes 

Tm OC r*'/5-3/5 (3.29) 

for discs with general A/(r) (see equation (|3.6p ). Tm oc r~^'* 
for the MMSN disc models (s = -3/2) while Tm oc r"^ * 
for the disc models with s = — 1. These steep profiles obvi- 
ously indicate that planets migrate outwards. For the region 
around the ice lines, R = RqT'"^ with Rq — 2 x 10^®. Conse- 
quently, equation (|3.14l) becomes 



Tm oc r 



s/5-3/20 



(3.30) 

for discs with M{r). Tm oc r""'^ for the MMSN disc mod- 
els (s — —3/2) while Tm oc r"" '* for the disc models with 
s = — 1. Under these shallow profiles, the entropy- related 
horseshoe drag cannot exceed the Linblad torque, and hence 
planets migrate inwards. Thus, a planet that migrates in- 
wards beyond the ice line will be halted near the ice line 
due to this opacity transition. 

3.3 Comparison 

We compare the heat transition barriers to the ice line barri- 
ers. We especially focus on the relative location of each bar- 
rier. As discussed above. Fig. |3]shows the relative location of 
these two barriers (see the black and gray lines). The ice line 
barrier (equation p. 27^ 1 is located inside of the heat tran- 
sition barrier (equation (|3.26|l ') for discs with a wide range 
of Q. However, we again emphasise that the above analyses 
for these two barriers can be valid only for the disc region 
with (10~^ 1^ 01 1^ 10~^). In the region with a low value of 
a such as dead zones [ct ~ 10~^), any corotation torque will 
be saturated (i.e. zero), so that these two barriers are never 
activated there. We also note that for such low values of a, 
ice lines can become barriers as a consequence of the Lind- 
blad torque rather than the corotation torque. In addition, 
if the adiabatic approximation breaks down, which can arise 
for the late stage of disc evolution or the outer part of discs, 
these two mechanisms cannot be effective. We now examine 
the mechanisms of barriers activated by Lindblad torque. 



4 BARRIERS ARISING FROM LINBLAD 
TORQUES 

Let us now examine regions wherein any torques aris- 
ing from corotation resonances are negligible, so that the 
direction of migration is determined only by the Lind- 
blad torque. Such regions occur for discs with a shallow 
temperature pr ofile or a low value of a in dead zones 
(~ 10~^ HP10:|Paardekooper et al.ll201ll ; iMasset fc Casolil 
I2OIOI : iHasegawa fc Pudritzll201ll ). 



4.1 Lindblad torque 

4.1.1 Basic equations 

We adopt analytical formulae o f the L indblad torques de- 
scribed in iHasegawa fc Pudrit j l|201ll . see their Appendix 
A for the complete expression). I n the f ormul ae, the stan- 
dard Lindblad torque derived bv IWardI l|l997h is modified, 
so that the vertical thickness of discs is taken into account 
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Figure 3. The location of the water-ice line as a function of the 
turbulence parameter, a. For every panel, the solid and dashed, 
black lines denote the radius derived from viscous heating using 
So and 0.1 X Eq, respectively (see equation 113.2711 and Table (Sjl. 
The dotted lines are from stellar irradiation (see equation |I3.28| |V 
For comparison purposes, the heat transition radius is denoted 
by the two gray lines (see equation 13.2611 . For all panels, viscous 
heating defines the location of the water-ice lines. In addition, the 



without rigorously solving the 3D Euler equations. Thus, we 
treat the Lindblad torque density d?T^ /dzdr of a layer at 
the height z from the mid-plane which is written as 



(4.1) 



where e = +(— ) for the outer(inner) resonances, the 
wavenumber m is treated as a continuous variable; 



m(r, z) 



(4.2) 



f2p = ^ GM, I rj^ is the angular velocity of a planet with 
r = Tp, the angular velocity f2 is 



n"(r,2) = n|f,41 + ^ 



2\ -3/2 



rp dr ' 



the epicyclic frequency k is 



^ = mCa/vK, and the forcing function tp is 



(4.3) 



(4.4) 



2m ^ K ^ ^ 



where 



A(ar,C) = "1 



Ko{A) 



(4.5) 



(4.6) 



Or — r/rp is the resonant position, and C, = zjr (also 
see T able [H. We have adopted the standard approximation 
for V> (|Goldreich fc Tremainjl980l : |jang -Condell fc Sasselovl 
l2005h . In addition, we have added a factor f2/K in the second 
term of following MG04. 

Finally, the total torque exerted by the planets on discs 
is found by summing over all layers; 



dz 



■ d (dV^ 



(4.7) 



Note that the integration of equation (|4.1[) in terms of z 
with C = leads to the famous an alytical form ulae of the 
Lindblad torque density in 2D discs (|Wardlll997l ). Also, note 
that the tidal torques shown above are exerted by the planets 
on discs. Thus, the positive sign of e in equation (|4.1|l means 
that the planets exert the tidal torques on discs, resulting 
in the loss of angular momentum from the planets and vice 
versa. 

The validity of the above approach f or mimicking the re - 
duction of the ti dal torques in 3D discs (iTanaka et aLll2002l ) 
is confirmed bv lLubow fc Ogilvid (|l998 l) who showed that, 
in thermally stratified discs, 2D modes carry off more than 
95 per cent of angular momentum and vertical modes are not 
dominant for transferring an gular momentum between plan- 
ets and their disc (also see WardI Il988l ; lArtvmowicj Il993l ; 
Ijang-Condell fc Sasselovl lioOsF 
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4-1.2 Basic assumptions 

We assume discs to be geometrically thin (z = 0) and Kep- 
lerian. Under this assumption, we have 

E 

P - ^ 



(4.8) 



.3/2 



Kep •> 



(i + -yrT? 

V m 



Thus, the torque density becomes 



d fdr\ 



£2/1 GAL 



(1 + 4^2) H 



(4.9) 
(4.10) 

(4.11) 



where 
and 



i i + ^W.(A) + (^ + 2yrT? 



ifo(A) 



A : 



la,- - 11 



(4.12) 



(4.13) 



Thus, the behavior of the torque density depends on 
Yia^ijP /H. In the following subsections, we expand every 
quantity (Qr,S,_ff, and ^) as a series in l/m(^ 1), by as- 
suming simple power-law disc structures in § 14.21 and non- 
power-law disc structures in § 14.31 and 14.41 This expansion 
is assure d, since the torque takes the maximum value at 
10 (|Wardlll997l ). and consequently, allows us to derive 
simple relations governing the direction of migration. Thus, 
we perform local analyses as done in § 2. 



4.2 Simple power-law structures 

We derive an analytical relation for the direction of plan- 
etary migration. Here, we assume simple power-law struc- 
tures for the background density and temperature in discs, 
that is, E oc and T oc r'. As a result, the relation depends 
only on the exponent of E and T. 

We can further simplify the mathematics by taking the 
limit ^ — )■ (since the effect of gas pressure becomes negligi- 
ble in the Keplerian discs; see Appendix|B]). Equation (I4.11|) 
then becomes, to first order in 1/m, 



dz \ dr 
where 



2e f t 1 5 
1+3^^^-2 + 2 + 4 



ro(rj,) = 2/i2GAf.§^mVo 



Thus, the total torque is 



dzTo{rp) 



3m 



(4.14) 
(4.15) 



3^-ro(r,)iJ,^s-i + Z 



(4.16) 



where the integration range for z is approximated to extend 
over the height of the photosphere Hp ~ 2Hp at the loca- 
tion of the planet since the density above the photosphere 
is about two orders of magnitude less than that at the mid- 
plane. 



Consequently, the sign of the net torque depends on the 
sign expression; 

t 7' 



sgn s 



2 4 



(4.17) 



The torque becomes positive (inward migration) when s — 
t/2 -I- 7/4 > and negative (outward migration) when s — 
t/2 + 7/4 < 0. When the MMSN disc models (s = -3/2) 
are adopted, t> 1/2 is needed for outward migration while, 
when disc models have s = — l,f>3/2is required. In 
general, T is a decreasing function of r. Therefore, some 
specific features in discs are required for the Lindblad torque 
to excite barriers. In the following subsections, we focus on 
dead zones and ice lines. 



4.3 Dead zone barriers 

Dead zones are known to play very i nteresting roles in plan- 
etary formation and migration (e.g. iMatsumura fc Pudrit j 
,2006) . They can excite two barriers: density barriers 
(MPT07) and thermal barriers (HPIO). These two barri- 
ers are interesting products of turbulent transitions. In this 
subsection, we derive the analytical relation which predicts 
when the direction of migration is reversed from inwards to 
outwards for these cases. 



4-3.1 Thermal harriers 

Thermal barriers were uncovered by HPIO who first investi- 
gated numerically the effects of dead zones on the tempera- 
ture structure of discs (also see[Hasegawa & Pudritz 2010t|). 
In their paper, the temperature structure of discs with dust 
settling and a dead zone was simulated by solving the full 
wavelength dependent, radiative transfer equation by means 
of a Monte Carlo method. Stellar irradiation is assumed to 
be the main heat source. HPIO demonstrated that planets 
migrate outward in a region where a positive temperature 
gradient is established (see their figs 1 and 3). This temper- 
ature gradient is a result of the back heating of the dead 
zone by a thermally hot dusty wall, and is well represented 
by T oc r* with t' > 3/2. The hot dusty wall is produced by 
the enhanced dust settling in the dead zone and the resul- 
tant enhanced absorption of stellar irradiation by the wall. 
Since they adopted E oc r'' with s — —1, the profile of the 
positive temperature gradient is identical to the profile re- 
quired by equation (|4.17|) in order for planets to migrate 
outwards. Thus, our relation predicts the results of the de- 
tail numerical simulations very well. 



4.3.2 Density jumps 

Density barriers were found by MPT07 who undertook a 
pioneering study on the effects of dead zones on planetary 
migration. They are created by the formation of a density 
jump at the outer edge of dead zones which arises as a con- 
se quence of time-de pendent, viscous evolution of discs (also 
see lMatsumura et al . 2009, for a more complete study). This 
density jump develops because of the huge difference in the 
strength of turbulence (a) between the active and dead zones 
of the disc. 

Here, we analytically derive an analogous relation for 
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density barriers. We adopt equation (|3.H) as the background 
surface density. In tiiis case, Fg is the density difference be- 
tween the active and dead zones, rtrans = r^dge is a orbital 
radius of the outer boundary of the dead zone, and oi = cH 
is the width of the transition (also see Table [l]). Except for 
this, we adopt the same approximation as above. As a result, 
the relation depends on the structure of the dead zones as 
well as the exponent of S and T. Since planets are slowed 
down or stopped near r^dge, we can approximate r^dge ~ rp. 

We note that a slightly more accurate shape of the 
boundary in equation (|3.1|l can be expressed in terms of 
error functions. However, the difference between the error 
and tanh functions i s less than 1 per cent with the proper 
choice of coefficients l|Bassetll 19981 ). Furthermore, tanh func- 
tions are fully analytical. For these reasons, we adopt tanh 
functions. 

Adopting the same approximation as above, we find 

+ hp- 



Table 5. Summary of parameter study 



1 + ^ 



1 — tanh 



2c cosh''(e/c) 



hp 



(4.19) 



(4.18) 

Note that we expanded E in terms of 1/m and finally ex- 
pressed it in terms of hp assuming 

2 
3m 

This is reasonable because the resonant positions are pushed 
away from th e planet by a distance ~ H(= rh) due to the 
gas pressure ()Artvmowiczlll993l : IWardlll997l . also see equa- 
tion groj). 

Inserting the above surface density into equation (14. lip , 
the sign of the net torque consequently becomes 



sgn 



-{Fg - l)tanh 



+ (Fg + l)hp s - 



t 



(4.20) 



Again, we can check the validity of equation (I4.20|) by 
comparing it with the numerical simulations. In this case, 
we compare the results of MPT07. Changing parameters Fg 
and c in equation (|3.1|) . they investigated what values of 
Fg and u result in outward migration (see their Appendix 
B). Table [5] summarises their experiments with our predic- 
tions (given in the brackets) . One observes immediately that 
equation (|4.2U|I well explains the results of their numerical 
simulations. One might wonder if the corotation torque is 
significant around the density jump, but this is not the case. 
MPT07 clarified that planets migrate outwards due to the 
reversal of the balance of Lindblad torques, not corotation 
torques. This arises because the Lindblad resonant positions 
are located further away from planets than the corotation 
ones. Therefore, migrating planets first encounter the inner 
Lindblad torque at the density jump and this is sufficient to 
reverse the direction of migration. 

In summary, both our simple analytical relations (equa- 
tions (|4.17|) and (|4.20p ) well reproduce the results of the de- 
tailed numerical simulations of dead zones as planet traps. 
This indicates that our assumptions and treatments are rea- 
sonable to capture the physics arising in the more compli- 
cated numerical simulations. 

4.4 Ice line barriers 



We finally examine ice line barriers. As discussed in § 13.2.41 
the disc radius of the ice lines is determined by viscous heat- 



Fg=5 


Fg=10 


Fg=100 


c = 4 in (0.07) 


in (-0.28) 


out (-6.6) 


c = 2 out (-0.8) 


out (-2.2) 


out (-28) 


c = 1 out (-2) 


out (-4.9) 


out (-57) 



/ip=0.35, s=-3/2, and t=-l/2. in(out) means the 
inward(outward) migration calculated by MPT07. The number 
in the brackets is calculated from equation II4.2UI1 . 



ing for a wide range of a (see Fig. [3]). In this subsection, we 
discuss how ice lines establish layered structures, and de- 
rive analytical relations for the resultant barrier. Thus, we 
investigate ice lines as an example of opacity and turbu- 
lent transitions. Again, we specify a ice line due to water, 
although all our analyses are applicable to ice lines of any 
material. The discussion is presented in § 15.31 



4-4- 1 Layered structures at ice lines 

As mentioned before, the number density of icy grains sud- 
denly increases at the ice lines. If discs are considered to 
be turbulent due to the MRI, this sudden increment of dust 
can result in layered structures: the MRI active , surface layer 
and the MRI dead, inner layer (|Gammiell 19961 ). This can be 
understood as follows. At the ice lines, the number density 
of free electrons in the disc also suddenly drops because they 
are absorbed by such dust grains. Since free electrons are the 
main contrib utor of coupling with magnetic fields thread- 
ing the disc (|Sano et al.l |200(]| ). the surface density of the 
MRI active region is strongly diminished, and consequently 
the surface density of the MRI dead region is enhanced. As 
a result, a density bump appears which acts as a barrier 
(see Fig. [4|. Thus, the mean value of a around the ice lines 
becomes a ~ an (see equation (14.21^ ). and hence, the ice 
lines can be regarded as a self-regulated, localised dead zone. 
In summary, the subsequent process initiated by icy grains 
can reduce the value of a at the ice lines. The disc radius 
at which it appears is determined by viscous heating (see 
equation p.27|l ). 

We recall from § 13.2.41 that ice lines become a barrier 
due to the corotation torque (also see IL08). However, any 
corotation torque may saturate at the ice lines because the 
turbulence dies out there - a conclusion which is valid for 
any disc in which turbulence is excited by the MRI. If the 
MRI is the most general source of disc turbulence, we expect 
that the ice lines may be natural barriers arising from the 
Lindblad torques rather than corotation torques. 



4-4-S Disc models 

We adopt a parameterised treatment of dead zones in a disc 
described by IL08. In the region with layered structures, the 
effective a can be generally written as 



Eaqa + (E — Ea)qo 



(4.21) 



where Ea is the surface density of the active layer, and aA 
and au are the stren gth of turbulence in the active and dead 
layers, respectively (|Kretke &: LinI l2007l : iMatsumura et all 
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I2009I . also see Table [T]). Assuming stationary accretion discs 
(see equation p.Gp ') with equation (|4.21[) . the surface density 
is given as 



M 



S = ^^-E.^A^, (4.22) 

where Ea cannot exceed E. Thus, equation (|4.22ll is useful 
for representing the surface density within dead zones where 
Ea < E. Also, this equation implies that the structure of E 
in the dead zones strongly depends on T,a- However, the 
structure of Ea is not well constrained by theory, since it 
is very sensitive to the distribution of dust grains as well 
as th e chemical models (|Sano et al. I I2OOOI : lllgner fc Neisoiil 
l2006f). Therefore , we a dopt simple prescriptions proposed 
bv lKretke fc LinI l|200'it and IL08 as follows. For the active 
region where Ea ~ E, the surface density simply becomes 

E^^^, (4.23) 

if there is no ice line. The effects of the ice lines on E are 
described below. 

4.4-3 Density jumps without ice lines 

At first, we examine density barriers produced by dead zones 
(see equation (I4.22|) ) to compare the our previous analysis 
done in § 14.31 For the surface density of the active layer, we 
assume 



Ea — E AO /ice — 



(4.24) 



following iKretke fc LinI (|2007l ). We note that Eao and /ice 
are very sensitive to the dust distribution (also see Table [T]). 
However, we set Eao constant and /ice = 1 here. A more 
detailed analysis is presented below where the effects of ice 
lines on /^ce are included. In discs with layered structures, 
it is generally considered that sa ^ 0. 

Using the above equation, equation (I4.22p is approxi- 
mately written as 



M 



, CIA 



1 — ehr, 



3 

*+2 



(4.25) 



-EA(rp)- 



O-D 



Q.D 



(1 + esAhp) 



where we performed Taylor expansion of E in terms of 1/m, 
and finally expressed it by hp, assuming hp ~ 2/3m. Thus, 
the direction of migration is determined by 



sgn 



M 



SnanHpflp 



-h„ 



3 1 

2^+4 



(4.26) 



-EA(rp 



^ OA — ao 



hn 



SA 



t 7 
2 + 4 



Since the outer edge of dead zones where a density jump 
develops is established around the region with Ea ~ E/2 , 
we can assume 



EA(r 



^ QA — «D 
CtD 



M 



OLA — OLD 



(4.27) 



Z-KCLdH'^^P OLA + OLD 

We emphasise that equation (|4.27[) gives the position of plan- 
ets halted by the density barrier which defines the outer edge 



of the dead zones (vp 



' edge ) 



In addition, this equation is 



another expression that the necessary condition (Ea < E) 



is safely satisfied. As a result, the sign of the net torque is 
written as 



sgn 



t 7 

Sa \ — 

2 4 



CtA + CtD 

2a.D 



SA + t + - 



(4.28) 



Compared with equation (|4.17() . the layered structure with 
a density jump gives an additional term which controls the 
direction of migration (since a a old). Thus, planets mi- 
grate inwards when sgn(sA -l-i-l-3/2) < while they migrate 
outwards when sa 4- i 4- 3/2 > 0. In standard disc models, 
SA > and t ~ —1/2. Our relation therefore dictates out- 
ward migration, which is also consistent with our previous 
analysis done in § 14.3.21 



4.4-4 Density jumps with ice lines 

We incorporate the effects of ice lines on the surface density 
of the active layer, following IL08 (see Fig. 14)1. In this case, 
we need to consider two separate cases, depending on the 
ratio of ru to Vedge, where ru is the disc radius of an ice 
line. For the case that ru/redge < 1, we can use equations 
(|4.22|l and (|4.24p as the surface density while we need to 
modify equation (|4.23|l for the case that ru/redge > 1. 

We first examine the case that ru < redge, that is, an 
ice line is located inside of a dead zone. In this case, /ice in 
equation (14.241) is generally given as 



fice — fd 



l + Z^lfl + tanhf^ 
2 V \ w 



0) 



(4.29) 



where fd and fai parameterise the effects of an ice line at 
which the s urface density of the active layer suddenly drops. 
As noted bv lKretke fc LinI (|2007l 'l. the surface density distor- 
tion induced by ice lines produces a radial, positive pressure 
gradient there, and consequently dust can be trapped there. 
Otherwise, it migrates inwards due to the so-c alled head 
winds arising from the gas motion (|Weidenschillin g _1977i). 
Since the ma:ximum migration speed of dust is ~ 1/a times 
faster than the viscous evolution of gas, dust is quickly piled 
up at the ice lines. This effect can be included in fd which 
is generally written as 



fd = 



1 + fd2 Cxp - 



ru 



(4.30) 



We note that inclusion of the dust effects may upgrade the 
ice line barriers to planet traps although we keep calling 
them barriers. 

Fig. |4] shows the typical E (as controlled by equations 
H4.29|l and (|4.30p ) and the role of the parameters fdi and 

//2 in /ice . 

Consequently, we can approximate the surface density 
(see equations (fi^ . (g^H), (gjS}, and (gjO])) as 



M 



SivH^QpaD 

^A^r'p, fice - 

where 

D = AB + ehp 

B [AsA + e 



ehp t + 



(4.31) 



ao 



t + 3 fdi 1 
2c 2 cosh2(l/c) 



A{B-l 



(4.32) 

t-f 3\ 
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e tanh ( - 



and 



B = 1 + fd2 exp — - 



(4.33) 



(4.34) 



We emphasise that the location of the ice line barrier which 
can halt planets migrating inward (rp « m) is determined 
by viscous heating (see equation (|3.27|l ') while the location 
of the density barrier produced by dead zones (rp « r^dge) 
is determined by equation (|4.27|l . 

As a result, the sign of the net torque is given as 



sgn 



M 



3 1 

'2^+4 



^A{rp, fi, 



(4.35) 



where 
E = A 



(4.36) 



+ hpAi 



-{B-r 



t + 3 



t 7 

-2 + 4 



The above equation is identical to equation (|4.27l) if A\ and 
-B — ^ 1, and A2 — > (equivalently, fdi and fd2 — > 0). Be- 
fore fully examining equation (|4.35p which looks very com- 
plicated, it is very useful to understand a simplified case. 
Taking the leading terms, outward migration arises when 



^AiVp, fice = 1) > 



M 



3Tv{aA — an) rpCs{rp) 



(4.37) 



The relation shows that the ice line barrier is readily active 
for discs with low accretion rate, high disc temperatures, 
high surface density of the active region, and large disc radii 
of the ice lines {rp ~ ru). In addition, this relation is very 
useful for constraining the effectiveness of the ice line barrier, 
which is discussed more in the next subsection. 

We can readily solve the equation (|4.35|l . since it is 
quadratic in t. Fig. [5]shows ru as a function of t. The black 
lines denote the negative solution while the gray, thick lines 
are for the positive solution. The regions where planets mi- 
grate outwards are encompassed by one of the solutions, 
X— and, axes (labeled by Outward). In this figure, we 
set fdi = 6, fd2 =2, SA = 3, and Eao/Sq = 0.01 at ro = 1 
au following IL08 (also see Table [3l l. We consider a water-ice 
line, so that T tt oiru) = 170 K (Ijang-Condell fc Sasselovl 
200J). We use c — I, because an y sharp transition in d iscs is 
smoothed out by disc viscosity (|Yang fc Menoull2010h . For 
stellar parameters, we adopt the values of CTTSs (see Ta- 
ble [3)1 . This disc setup is our fiducial model (see the solid 
lines). Based on equation (|4.37p . the choice of stellar param- 
eters (i?» and Tt) does not matter if is parameterised (see 
equations (|3.24|) and (|3.27|l '). Although the complete treat- 
ment in which ru depends on stellar parameters is presented 
in § [S] we simply label ru by asterisks for comparison pur- 
poses here (see equation (|6.1[) ). 

In Fig. O we explore parameter space by changing M, 
T,Ao, and sa (on the top, middle, and bottom panels). For 
t <0 which is general for any disc, higher M and sa, and 
lower Eao reduce the outward migration region. This is con- 
sistent with the above argument (see equation 1)4. 37p '). The 
situation is the opposite for t >0 which can happen due 
to hot dusty walls. If the dependence of ru on Al is taken 




1 00.0 



Figure 4. The typical structures of S and controlled by fice- 
The upper solid line denotes the case that an ice line is located 
within a dead zone while the lower solid line is for an ice line 
beyond the dead zone. The distribution of S 4 is denoted by the 
dotted line. The definitions of f^i and /;j2 in fice and of in S 
are labeled in the figure. 



into account (see the top panel), higher accretion rates do 
not necessarily reduce the region where the water-ice line 
barrier is active. This arises because higher accretion rates 
result in larger water-ice line radii which are preferred for 
outward migration (see equation (|4.37p i. Thus, the condi- 
tions for outward migration are controlled by the complex 
dependence on M. Another interesting feature is shown in 
the bottom panel. Around ru ~ 30 au with t > 0, both 
negative and positive solutions merge with each other (see 
the black and gray dotted lines on the bottom panel). This 
happens because there is no solution to equation (|4.35p . 

Now, we examine the case that ru/redge > 1, that is, an 
ice line is located beyond the dead zone in the active region 
of the disc. In this case, we adopt the following functional 
form for E; 



M 



3TvcsHaA 



1 + gd exp 



■ ru 



(4.38) 



Compared with equation (|4.23|l , the ice line produces a den- 
sity bump at r = ru which is a consequence of a localised 
layered structure. This function well represents the results 
of IL08 with a proper choice of gd (see the lower solid line 
in Fig. [4]). Here, we only examine the power-law index of E 
rather than expanding it in terms of h. This is because one 
of our approximations, rp ~ ru, limits the applicability to 
Gaussian-functions. In fact, planets migrating towards the 
density bump can be halted around rp = ru+ A r where a 
negative steep surface density profile plays the critical role. 
However, the assumption that rp ~ ru washes out this effect. 
This also indicates that equation 1)4. 35|) may underestimate 
the required condition. 

The power-law index of E is given as 



dluE _ 3 
^ ^ dlnr ^ ^ ^ 2 ^ 



'-i gd exp 



1 + gd exp - 



(4.39) 

We recall that s becomes equivalent to s for pure power-law 
discs (see equation (|3.2|) ). Assuming planets to be halted 
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Table 6. The condition for the ice Une barrier to be active 
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Figure 5. The location of a water-ice line as a function of the 
temperature exponent t. The black lines denote the negative so- 
lution to equation 114.351 1 while the gray, thick lines are for the 
positive solution. The regions where planets migrate outwards are 
surrounded by one of the solutions and x— and jy— axes (labeled 
by Outward). For the top, middle, and bottom panels, M, T,aO} 
and sa are varied, respectively. Otherwise, we adopt the values 
of the fiducial model. For the case of t < with higher M and 
sa and lower S^o> outward migration region shrinks. For the 
case of t > the situation is the opposite. For comparison pur- 
poses, we label ru by asterisks (see equation 116. Il l V For the top 
panel, the value of ru increases with increasing M. This indicates 
that the conditions for outward migration are determined by the 
comolex deoendence on M (see eauation (14. 3711 '1. 



/In 2 where the Gaussian function takes 
half of the maximum value, the resultant value of s becomes 



l + 5d/2 



(4.40) 



As a result, we find that the direction of migration is con- 
trolled by 



sgn 



(4.41) 



where we have assumed any disc quantity to behave as 
power- law in a local region centered at r ~ rp (see §[4]2|. For 
discs with hp — 0.1, c = 1, and gd ~ 3, t > —6.5 is required 
for outward migration while for discs with hp = 0.1, c = 1, 
and gd = 6, t > —8.2 is needed. Thus, the direction of mi- 
gration is readily reversed by the density bump produced by 
the ice lines. Furthermore, larger density bump gd expands 
the outward migration region, as expected. 



4.5 Locations of the barriers 

We can now compare the dead zone and ice line barriers in 
order to examine the relative importance of each. We adopt 
the analysis done for the case of ru/redgc- We especially 
focus on the relative location of each barrier. Combining two 
equations (I4.27|) and (|4.37|l . we find the condition which is 
given as 



> h{redge) 



OA + ceo \ «7+*72+T 

OtA — CtD 



(4.42) 



where we set ro — Vedgc. If this condition is satisfied, the ice 
line and density barriers are both active while if not, only 
the density barrier is active. Thus, the separation between 
the ice line and dead zone barriers is required to be relatively 
small in order for the ice line barrier to be effective. Table [6] 
summarises the threshold values of ru/r^dge for various sa, 
t, and h{redge). As all three quantities increase, the threshold 
values of ru/r^dge increase. In other words, the parameter 
space where the ice line barrier is effective is larger for lower 
values of sa, t, and h{r^dge). More complete treatments in 
which the ice line radius is determined by viscous heating 
are presented in § |S] 



5 DISCUSSION 

Before we integrate our analyses into a coherent picture, 
we discuss other possible barriers induced by other physical 
processes which are neglected in this paper. Also, we discuss 
opacity transitions that are neglected in our unified picture. 
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In addition, we examine the possibility of molecules other 
than water to excite ice line barriers. Finally, we estimate the 
mass range in which planets are regarded as type I migrator. 

5.1 Other possible barriers 

We have so far assumed that stellar irradiation heats up 
only the dust in discs. However, photons emitted from stars, 
especially with very short wavelengths such as UV and X- 
rays can heat up the gas. This resul t s in the photoevapora - 
tion of discs (HoUenba ch et al.|[l993 : Ijohnstone et aLlliggsl ) 
- wh ich co nsiderably affects planet formation a nd migration 
fe.g. iMordasini et all [20091 ). iLvra et al.l (|2010l ) showed that 
photoevaporation can activate another barrier. This arises 
by the reduction of the surface density of gas due to photo- 
evaporation. Since viscous heating strongly depends on the 
surface density (see equation (|3.14|l '). this reduction results 
in shallower temperature profiles. Consequently, planets are 
halted at the region where the transition of the temperature 
slope occurs. 

In general, photoevaporation is considered to be im- 
portant for discs surrounded by nearby OB associations 
which provide a huge input of high energy photons. In ad- 
dition to photoevaporation, one expects that such massive 
stars can significantly affect the disc structure by changing 
the dust temperature. Recently, iGorti fc HollenbachI (|200£ ) 
have shown that, without nearby massive stars, even low- 
mass central stars (> O.SMq) can drive photoevaporation of 
their surrounding discs (r > 20 au). This arises from ener- 
getic stellar irradiation (far-UV, extreme- UV, and X-rays). 
We will address this in a future publication. 

In addition, we neglect the effects of the inner edge of 
dead zones where a positive radial gradient of the surface 
density can form. Such a profile can produce a barrier due to 
the vortensity-related corotation torque. However, the disc 
radius of the inner edge of dead zones can be fixed around 
r ~ 0.01 au, because the radius is controlled by the thermal 
ionisation temperature of gas. Thus, the barrier is unlikely 
to play an important role in the diversity of the detected 
exoplanetary systems. For the reason, we disregard it in this 
paper. 

5.2 Neglect of opacity transitions 

We have intensively investigated ice lines as an opacity tran- 
sitions in § 13.2.41 while our approach is applicable to any 
other opacity transitions. In general, protoplanetary discs 
have several opacity transitions. MG04 showed that any 
opacity transition (including ice lines) affects the disc surface 
density and temperature in a similar fashion. Therefore, it 
is naturally expected that a single disc can have several bar- 
riers that are all excited by opacity transitions. Nonetheless, 
we will neglect all opacity transitions in [|6]for the following 
reasons. 

Except for those produced by ice lines, all opacity tran- 
sitions are a consequence of the high disc temperatures and 
subsequent destruction of opacity sources. As a result, their 
disc radii distribute well inside of 1 au (see fig. 1 of MG04). 
MG04 found that these opacity transitions make the mi- 
gration time significantly longer due to the Lindblad torque 
alone. (Equivalently, the net torque exerted on a planet be- 
comes a very small (in magnitude), negative value.) This 



indicates that, in order to actually stop the planet's inward 
migration, a small amount of outward-directed corotation 
torque would be needed. However, any corotation torque at 
the opacity transitions is likely to vanish due to the presence 
of dead zones. As discussed above, protoplanetary discs may 
have dead zones that are 1- 10 au in size from their central 
stars, depending on their column density. One of the key fea- 
tures of the dead zones is a low level of turbulence, which re- 
sults in the saturation (ineffectiveness) of corotation torques 
there (see §|4ll. Thus, opacity transitions within dead zones 
cannot be planet traps and therefore we disregard opacity 
transitions other than those created by ice lines. For the ice 
lines, it is more likely that the process initiated by the incre- 
ment of icy dust grains reduces a level of turbulence there 
and hence ice lines are regarded as localised dead zones (see 
I4.4.1|l . Thus, we include ice lines as an opacity and turbu- 
lent transition rather than an opacity transition in § [6] In 
summary, any opacity transition is neglected below. 



5.3 Ice lines of other molecules 

We have focused on a water-ice line in the above dis- 
cussions, although our analyses apply to ice lines of any 
molecules in discs. It is well known that water is not the 
most abundant chemica l species in protoplanetary discs (e.g. 
lAikawa fc Herbstlll999l ) In general CO is the most abun- 
dant molecule in discs. iDavid l|2007l ) showed that the total 
amount of CO-ice is much larger than water-ice in irradiated 
accretion discs. Furthermore, differen t spec ies have differ- 
ent condensation temperatures. iDavij ^2005) demonstrated 
that the location of ice lines strongly depends on molecules 
in consideration. Thus, it is very interesting to investigate 
which molecules can produce barriers in discs. 

In order to proceed, we adopt the analyses done in § 14.41 
We especially focus on the minimum abundance of molecules 
which is required for their ice lines to be a barrier. We 
present the detail analysis in Appendix [C] and summarise 
our findings here. We find that whether or not ice lines of 
other molecules act as barriers strongly depends on the disc 
parameters such as sa and t. Adopting reasonable values 
for them, we find that, for molecules that condense within 
a dead zone, only specific, abundant species such as CO 
have the possibility to work as a barrier. If molecules freeze 
out beyond the dead zone, then even a tiny fraction of the 
molecules can act as a barrier. In summary, other molecular 
species can produce ice lines but these will form near to the 
position of the outer edge of the dead zone (r^dge)- We leave 
more comprehensive models for ice lines to a future publica- 
tion. In the following discussions, we again assume water-ice 
to be abundant enough for creating a barrier. 



5.4 Mass limits for type I migration 

We discuss the mass limits that our analyses apply to. Type 
I migration is only applicable to low mass protoplanets. 
Massive protoplanets can open up a gap in their discs and 
undergo type II migration. This mode of migration is con- 
trolled by viscous evolution of discs and is much slower than 
the type I migration. The critical mass, also known as the 
gap-opening mass, is well discussed in the literature (e.g. 
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iMatsumura fc Pudrit dlioO^ '). and given by 



^3Chl,C^Jmahl^, (5.1) 



where we have introduced a coefficient C — 3 into the orig- 
inal condition by taking into account the reduction of the 
tidal torque due to the disc thickness l|Hasegawa fc Pudrit j 
The left term in the brackets of equation (15. 1[) is 
derived by a Hill radius analysis while the right term for 
viscous ones. Table [7] summarises M^^max for stars with var- 
ious masses at rp = 1 au. The stellar and disc parameters are 
given in Table [S] For discs with a — 10~^, the critical mass 
is established by the Hill radius analysis, and has the same 
order of Neptune to Saturn masses. On the other hand, the 
criterion derived for viscous discs provides the gap-opening 
mass for discs with a = 10~^, and Mp^max is on the order 
of Earth masses. 

What is the minimum mass for protoplanets to un- 
dergo type I migration? It is interesting that the dynamics 
of low mass bodies, from dust grains to plan etesimals, can 
drastically differ f rom that of protoplanets (^ Adachi et al.l 
ll976l : lBai fc Stone 2010: Nelson & Grcsscl 201(|). This arises 
mainly from gas drag. One of the most famous conse- 
quences of the gas drag is the rapid inward migration of 
meter-sized particles discussed in § 14.4.41 The size depen- 
dency of gas drag affects the formation of planetesimals (see 
IChiang fc Youdin 2010, for a recent review). However, we 
simply focus on the dynamics of planetesimals, since we are 
interested in minimum mass of type I migrator. In laminar 
discs, the rate of cha nge of semi-major axis induced by gas 
drag is given by (see lAdachi et aLlll976l . for derivation) 



drp 
~dt 



drag ~ 



M, 



1/3 



(5.2) 



where Kdrag is a function of the eccentricity, the inclination, 
and the material density of a plane tesimal, and th e radia l 
pressure gradient of discs. Following iKokubo fc Idal (|2000l ). 



we find Kdra 



2.7 X 10 in cgs units. As a result, the 



timescale of gas drag is written as 



^drag 



dt 



M, 



1/3, 



(5.3) 



where we have assumed pp ~ T,p/Hp. On the other hand, 
the migration timescale is generally written as 



Alp Vp ^p 
2f ' 



where the tidal torque is scaled by 



Mp 



^pTp^p 



(5.4) 



(5.5) 



and Kmia is an order of 1 t o 10, depending on disc mod- 
els jPaardekooper et al.ll2O10l ). Equating equation (|5.3|) with 
equation H5.4[) . the critical mass is given as 



Mp 



KdraghpM^ 
'^Kmig fp 



3/4 



(5.6) 



Table [7] summarises the minimum mass of type I migrator 



for stars with various masses. We set K„ 



5. The values 



of Mp.min are at least three orders of magnitude smaller 
than Earth masses. Thus, the dynamics of solids in laminar. 



gaseous discs derives from type I migration over four to six 
orders of magnitude in mass. 

In turbulent discs, the situation becomes more com- 
plicated, because planetesimals are affected by gas drag 
as we ll as stochastic torque a rising from density fiuctua- 
tions. iNelson fc Gressell (|2010l ') undertook numerical simu- 
lations of the dynamics of planetesimals in MHD turbulent 
discs and found that the stochastic torque becomes dom- 
inant over the gas drag for planetesimals which are ~ 25 
meters in size. This implies that the critical mass Alp.min 
is involved with the stochastic torqu e rather than t he ga s 
drag. Adopting the torq ue formula of iTanaka et al.l (|2002h . 
iNelson fc Gressell (|2010l i estimated that Mp,™i„ ~ 0.1 IM® 
in discs with a ~ 10~^. 

Our value for Mp^min, is estimated for laminar discs and 
is several orders of magnitude smaller than that in turbulent 
discs. However, equilibrium states for planetesimals larger 
than 100 m are not achieved in their simulations. Thus, 
more effort is required for accurately estimating the min- 
imum mass of type I migrator. 

In summary, our analyses are likely applicable to proto- 
planets with the mass range from well below Earth to Saturn 
masses for either laminar or turbulent discs. It is important 
to emphasise that almost observed low-mass exoplanets are 
covered by our analyses. 



6 A UNIFIED PICTURE OF PLANET TRAPS 

We now integrate the analyses of § O U and [S] into a sin- 
gle comprehensive framework for understanding planetary 
system architectures. More specifically, we investigate under 
what conditions which barriers are effective, and identify the 
location of the active barriers in the entire disc. Thus, we 
call any barrier a planet trap, since it plays a role of halting 
migration as well as of affecting the subsequent formation of 
planetary systems. We summarise all possible planet traps 
and the necessary conditions for them in Table [S] We find 
that up to three planet traps can exit in one disc model 
while the minimum case is just a dead zone trap. 

Before examining the location and conditions in Table 
[S] we re-derive equations (|3.27|l and p.26p by simply sub- 
stituting Eo with M using equation (|3.6|l . Then equation 
(|3.27|l becomes 



In 

ro 



27 Kg IIgQ,Q 



.H^O 



M 



(ru) GAasBCffkB \3-k 



2/9 



(6.1) 



where ko = 2 x 10 , and equation (|3.26p becomes 



rht 
ro 



1 



TmO V R 



ro \^^'^ ( 27^0/19^0 



M 



1/31 

(6.2) 



where iio ~ 2 x 10~*. Also, equation (|4.27[1 is rewritten as 
M 



•edge 

ro 



i+t+3/2 



(6.3) 



where we assume that sa, t, and Eao are all external pa- 
rameters. This rearrangement of equations is very useful for 
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Table 7. Mass limits of type I migrator at 1 au 



Herbig Ae/Be stars CTTSs M stars 



10-2) 63.8 50.7 18.6 

10-5) 3.0 1.9 0.6 

3.7 X 10-3 4.6 X 10-* 4.8 x 10-^ 



Table 8. Summary of all possible traps 





Locations 


Conditions 


Ice line traps 


equation (13.2711 


equation (|4.42| only for ru < r^^ge 


Dead zone traps 


eauation l|4.27p 


N/A 


Heat transition traps 


equation (13.2611 


rht > redge 



integrating our separate analyses done in § [3]and|4]into one 
picture. 

The important result is that the locations of all three 
traps are mainly controlled by M, but with different power- 
law dependencies. This suggests that some traps can merge 
with each other due to time-dependent, viscous evolution 
of discs. This is clearly shown in Fig. [G] In this figure, the 
radius of all three traps is plotted as a function of M by set- 
ting (sA,i, Sao) = (3, — 1/2, Eo/lOO) using the value of Eq 
in Table [3] For stellar parameters, we adopt the values of 
Table |3l We calculate Ho, assuming stellar irradiation to be 
dominant. As discussed in the above two sections, the heat 
transition traps can be active for our models with a — 10"^ 
and their disc radii are represented by equation (16. 2p (see the 
lower dotted lines). For comparison purposes, we also show 
rht for models with a — IQ-^ (see the upper dotted line). 
The ice line traps can exist in both active and dead regions. 
Thus, their disc radii are expressed by equation (|6.1[) with 
a = 10~^ and IQ-^ for the active and dead regions, respec- 
tively (see the lower and upper dashed lines). The solid line 
is for the location of the dead zone trap (see equation (|6.3p ). 
Taking into account the conditions for them (see Table [S]), 
the locations of all active traps are shown by the black lines. 



6.1 A picture of forming planetary systems 
around stars with various masses 

We present a picture of how planetary systems are estab- 
lished based on planet traps and discuss how host stars affect 
the size of planetary systems formed around them. 

6.1.1 General picture 

Fig. [7] schematically illustrates this. The disc radius of each 
planet trap comes from Fig.[6]for the case of CTTSs (also see 
Tableini). The structure of surface density is plotted based on 
our disc models discussed in 14.41 Since we do not model the 
time dependency of M that gives the realistic value of sur- 
face density, we simply adopt the values of Eo and 10^ x Eo 
in Table[3]at M = W'^MQ/year and M = IQ-^MQ/year, 
respectively. These choices do not affect our picture. 

At the early stage of disc evolution (M = 
10~^ MQ/year), the mass of protoplanets (denoted by the 
size of dots) which are captured in each trap is small, and 



they are widely separated. As the disc accretes onto the 
host star, the surface density of gas and the accretion rate 
decrease, and the position of the planet trap region with pro- 
toplanets captured moves inwards. Concurrently, the growth 
of these protoplanets proceeds. Since the rate at which each 
planet trap moves differs (see Table the separations be- 
tween these trapped planets shrink. Also, the solid surface 
density in each planet trap probably varies in accordance 
with the variation of the surface density of gas. Therefore, 
the growth rate at the dead zones and ice lines may be higher 
than that at the heat transition. At some point, the grav- 
itational interaction between trapped planets could exceed 
planetary migration and planet-planet scattering effects be- 
gin to dominate. We leave the later stage of evolution of 
planetary systems into our forthcoming paper. However, we 
can speculate how planetary systems evolve as follows. 

If the trapped (proto)planets distribute themselves 
closely enough, the most massive planet can surely survive 
the planet-planet interaction. On the other hand, the other 
two trapped planets could scatter inwards or outwards, or 
could be ejected from the system, depending on the mass 
of the most massive planet. Thus, our picture can in prin- 
ciple predict how planetary systems are formed, and what 
differentiates the formation of multiple planets versus single 
planet in a system. 

6.1.2 Dependence on stellar masses and accretion rates 

We discuss the effects of stellar masses and accretion rates 
on the scale of planetary systems based on our picture. Ta- 
ble [9] summarises the position of each planet traps at two 
different values of M for stars with various masses in Fig. 
[6] One immediately observes that the scale of orbital radii 
of planets captured in planet traps is a function of stellar 
masses and disc accretion rates. For Herbig Ae/Be stars, 
planetary systems can extend from the order of au to ~ 10'^ 
au. For CTTSs, the size of planetary systems shrinks to ~ 
10 au. For M stars, planetary systems may be well confined 
within a few au. We discuss more on the dependence on the 
stellar mass by comparing with the observations below. 

6.2 Implications for observations 

We discuss the validity of our picture by comparing with 
the observations. Fig. [6] shows that the location of the dead 
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Table 9. The position of planet traps 









Herbig Ae/Be stars 

redge (AU) 


ni (AU) 


r-M (AU) 


M = 
M = 


10" 
10" 


4(M0/year) 
'^(^o/year) 


3.9 
1.2 


155 
20 


830 
47 








CTTSs 

redge (AU) 


■ni (AU) 


TM (AU) 


M = 
M = 


10" 

10- 


6 (A/0 /year) 
« (A/0 /year) 


3.7 
1.2 


11.7 
1.5 


42.4 
2.4 








M stars 

redge (AU) 


ni (AU) 


Tht (AU) 


M = 
M = 


10- 
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« (A/0 /year) 
9 (A/0 /year) 


4.2 
2.3 


N/A 
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Figure 7. Our general picture of tile formation of planetary systems. The disc radius of each planet trap comes from Fig. \E\ for 
the case of CTTSs. The position of each planet trap is summarised in Table |9] The structure of surface density is based on our 
disc model discussed in § 14.41 The dots denotes trapped protoplanets, and their size represents their masses. At the early stage of 
disc evolution (on the left panel), the accretion rate and the surface density are higher, and the separation between protoplanets 
which are captured at each trap is larger. Planetary masses are smaller. At the later stage of disc evolution (on the right panel), 
the accretion rate and the surface density become lower, and the separation becomes smaller. The planetary masses are now much 
larger. We assume that solid density at the outer edge of dead zones and ice lines is larger than that at the heat transition, based on 
the higher surface density of gas. It well illustrates how the orbital movement of planet traps at different rates affects the formation 
of planetary systems. 



zone traps for all cases is on order of 1 au (see the CTTS 
case). This is in accord with the observations which show 
stati stically that the mo st gas giants are found around 1 au 
(e.g. IWright et alll2009l ). Therefore, our finding may imply 
that the dead zones play a significant role in the formation 
and location of these planets. 

For the case of Herbig Ae/Be stars (the top panel), 
there are three traps for high accretion rates. It is interesting 
that the range of disc radii for the heat transition and ice 
line traps covers the (large) range of orbital radii (~ 10^ — 
10^ au) of massive planets around massive stars such as 
A st ars which are obse rved by the direct imaging method 
fe.g. lMarois et aLll2008l ). For a lower accretion rate, the heat 
transition and ice line traps merge. However, the position of 
this merger is well separated from that of the dead zone 
trap. Thus, at least two planets may survive the planet- 
planet interaction. This multiplicity trend for massive stars 
is confirmed again by the direct imaging method (although 



the detected semi-major axes are considerably larger than 
our prediction.) 

For the CTTSs case (the middle panel) , three traps ex- 
ist at first, and then as the disc accretion rate drops, the 
ice line and heat transition traps merge into the dead zone 
traps. As mentioned above, the interaction of three active 
traps arises around 1 au. Based on our picture, planets in 
multiple systems are less massive than single systems, since 
more massive planets provide more destructive effects on 
planetary systems. This prediction well explains the observa- 
tion trend that (apparently) single planets sta tistically have 
large r masses than planets in multiple systems (|Wright et al.l 
I2OO9I ) . More recently, the Kepler mission also confirmed the 
same trend (|Latham et al.|[2011^ . Thus, the interaction of 
planets in their traps is essential for a systematic under- 
standing of how planetary systems are formed. 

For the case of M stars (the bottom panel), only the 
dead zone trap is active at first. The location of this 
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Figure 6. The location and evolution of all three planet traps. 



1-9 



-■-1 1 



The heat transition traps with a 



10- 



and a 



10- 



denoted by the lower and upper dotted lines, respectively (see 
equation H6.2|l l. The ice line traps with a = and a = 10~^ 

are denoted by the lower and upper dashed lines, respectively (see 
equation 1 I6.II 1V The dead zone traps are shown by the solid line 
(see equation ijG.Sp ). Taking into account their active conditions 
(see Table |8ll, all effective traps are denoted by the black lines. 
For the Herbig Ae/Be star case, all three traps are effective for 



trap can explain the observed Super-Earth around M stars 
(|Beauheu et al.|[2006l ). As the accretion rate decreases, the 
ice hne trap appears and then disappears due to the con- 
straint derived from equation H4.42p (see Table [8|. In this 
case, it is more sensitive to when and what mass of planets 
are captured at each trap. 



6.3 Parameter study 

We perform a parameter study on sa, t, and Eao (see Fig. [8]). 
We adopt the values of the CTTSs case. Fig.[8]shows all the 
possible cases. The locations of ra and r^t are the same on 
all panels (since they are independent of these parameters, 
see equations (|6.ip and (|6.2[l ). The position of r^dge is scaled 
by Eao while the dependence of M on r^dgs is controlled by 
the sum of sa and t (see equation (|6.3|l '). One immediately 
observes that the presence of multiple traps in a single disc is 
more likely than only a single planet trap over a wide range 
of M. Hence, the interaction of planets which are growing 
and lodged in their traps can be considered as essential. As 
discussed above, the growth rate of trapped protoplanets 
which controls when the planet-planet interaction initiates is 
crucial. We will address this issue in the forthcoming paper, 
by taking into account the growth of protoplanet as well as 
disc evolution. 



7 APPLICATION TO OUR SOLAR SYSTEM 

In this section, we apply our picture to the origin of the So- 
lar system - in partic u lar to the Nice rnodel a s developed by 
iTsiganis et all (|2005h : lMorbidelh et al.l (|2005l ): [Gomes et all 



In this model, the dynamics of outer planets is in- 
vestigated by simulating planetary migration in a debris 
disc and the subsequent planet-planet interaction. Thus, 
the model assumes that all gas in the disc is depleted and 
planetary migration is caused by interactions with planetes- 
imals. However, these N-body based simulations depend on 
their initial conditions which are not justified. iMore recently, 
[Mats umu ra et al. ( 201Q ) investigated planetary migration in 
dissipating gaseous discs and showed that the semi-major 
axes are established by planetary migration in gaseous discs. 
Once gas is depleted, other orbital quantities such as ec- 
centricities are determined by the planet-planet interaction. 
Therefore, it is crucial to examine the initial conditions of 
the Nice model in discs where gas is still present, in order 
to validate their results. 

Fig. |5] shows the most plausible case for the initial con- 
dition of the Nice model in which Jupiter is placed at 5.45 
an while Saturn at ~ 8.2 an. We find it by varying the 
parameters (sA,t,SAo)- In this figure, we show the two ex- 
pected traps: the ice line and dead zone traps. These two 
planet traps may be the building sites of the two gas gi- 
ants (Jupiter and Saturn). The formation of the icy giants 
(Uranu s and Neptune) may originate from scattering pro- 
cesses iTh ommes et alj [19991 ) . This idea is also supported 
by the Nice model. Thus, the core of the assumption in the 
Nice model is the initial semi-major axes of two gas giants 
(Jupiter and Saturn; see the horizontal solid lines). Fig. [9] 
indicates that these two traps deliver the cores of two gas 
giants to the initial semi-major axes for the Nice model to- 
wards the end of the gas disc's lifetime, when the gas ac- 



The Origin of Planetary System Architectures. I. 



21 






100.0 r 



3 10.0 r--. 



(s*.t.S4j = (3. 1.5/0) 




100.0 



10"° '0"' 10"° 10" 

Qccretlor rate {M^^^y year) 



10" 



(s..t.S,o) = (B. 1.5/0" 



Heat Trorsltion 
Ice Line 
Dead Zone 



10"° -Q-/ ^Q-= ^Q-i 

accretion rate (M^^^,/ ye 



100.0 r 



10" 



(s„t,Sj = (5/4,-l/?,10^) 



Heat Transition 
Ice Line 
Dead Zone 




1 0"^ 1 0"' ' 0"° 1 0"^ 1 0" 
accretion rate (M^.^^/ year) 



Figure 8. Parameter study on the location and evolution of three planet traps around CTTSs. The definitions of all lines are the 
same as Fig. \U\ All possible cases are shown. Our choices of parameters are shown in each panel. The presence of multiple planet 
traps in a single disc is more common. Therefore, the interaction of the planet traps is likely essential. 
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cretion rate have fallen to M — 5.5 x 10~^ Mq /year (see 
the vertical solid line). At this moment, the disc still have 
considerable amount of gas (~ 1 per cent of So), and hence 
the initial condition adopted in the Nice model may be too 
simplified. 

Another interesting point is that Jupiter hovers around 
the ice line trap (see Fig. Since the dust density there 
is significantly enhanced, the formation of the trapped core 
of Jupiter speeds up considerably. This may be a reason 
why the mass of Jupiter is larger than that of Saturn. Re- 
cently, the significant improvements of the Nice model were 
achieved by incorporating planetary migration driven by the 
gas discs (|Morbidellill20ld . references herein). These studies 
showed that a key is the mass ratio of Jupiter to Saturn. If 
Jupiter is more massive than Saturn, the initial conditions 
adopted in the original Nice model are readily achieved. This 
is possible if there is outward migration of these two gas gi- 
ants, which according to the model is consequence of captur- 
ing them into the 2/3 mean motion resonances (MMRs). If 
such a mass ratio of these two gas giants is not attained, both 
planets migrate inwards too rapidly to survive. Although 
these trends of migration happen without any planet trap, 
no one can currently explain such a mass ratio. Thus, our 
analyses are likely important for giving at least one expla- 
nation for deriving the initial conditions for the Nice model. 
We will address more comprehensive roles of planet traps in 
the Nice model in a future publication. 

As the separation of these two traps decreases with de- 
creasing accretion rates, the two trapped cores can induce 
scattering of any other cores formed around them. Thus, our 
approach is very useful to investigate the formation of plane- 
tary system more self-consistently. In the subsequent paper, 
we will include the growth of planets in viscously evolving 
discs. 



8 CONCLUSIONS 

We have systematically investigated how inhomogeneities 
arise in protoplanetary discs and how they create planet 
traps where the direction of planetary migration changes 
sign: from inwards to outwards. We analytically investigate 
the various mechanisms of planet traps activated by both 
corotation and Lindblad torques. For the traps related to 
the corotation torque, we have simply examined the disc 
structures that lead to outward migration. For the traps re- 
lated to the Lindblad torque, we have derived the relations 
which predict when the torque reversal occurs. 

Our analytical relations derived from the Lindblad 
torque succeed in reproducing the results of the previous 
numerical studies done for the dead zones: the density and 
thermal barriers (i; l4.3|) . This confirms that our assumptions 
and treatments are valid for capturing the detailed, complex 
physics occurring in the simulations. Furthermore, we have 
shown that the ice lines can serve as barriers due to the 
Lindblad rather than corotation torque. In discs with MRI 
induced turbulence, the ice lines create the layered struc- 
ture, so that the strength of turbulence (a) reduces to the 
value of the dead zones. Thus, the ice lines can be considered 
to be self-regulated, localised dead zones. We have also dis- 
cussed the possibility that a single disc can have several ice 
line barriers which become effective due to different species 
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Figure 9. Application to the Nice model. Two planet traps exist 
in a single disc (as Fig. O . The horizontal solid lines denote the 
initial condition of the semi-major axes for Jupiter (r = 5.45 au, 
see the lower line) and Saturn (r x 8.2 au, see the upper line), 
respectively in the Nice model. The core of Jupiter is captured 
at the ice line trap while that of Saturn is at the dead zone trap. 
At M = 5.5 X 10~^ Mq /year , the trapped cores reach the initial 
condition adopted in the Nice model. 



such as water-ice and CO-ice. A more detailed understand- 
ing of the structure of ice lines is required for addressing 
these features. Finally, we have estimated the mass range 
over which our analyses can apply. In either laminar or tur- 
bulent discs, our modeling can be useful for understanding 
almost observed exoplanets. 

We summarise our major findings below. 

(i) We have shown that the heat transition of protoplan- 
etary discs from viscosity to stellar irradiation dominated 
heating, activates a new barrier, which we call the heat 
transition barrier (see i; l3.2.3|) . This barrier arises due to the 
entropy-related horseshoe drag and is caused by the temper- 
ature transition there. At disc radii inside the heat transition 
radius, r^t, the temperature slope determined by viscous 
heating is very steep, resulting in outward migration. On 
the other hand, a shallower temperature profile controlled 
by stellar irradiation dominates at larger radii beyond r^t- 
Such a region drives inward migration. Thus, planets which 
migrate toward the location of the heat transition are cap- 
tured there. 

(ii) We have also shown that a single disc can have up to 
three planet traps: the dead zone, ice line, and heat transi- 
tion traps (see § |6} . This multiplicity of planet traps is likely 
quite common. We have derived fully analytical relations for 
their disc radii and their effective conditions. These relations 
imply that dead zone traps may be the most important in 
order to explain why so many gas giants are observed around 
1 au. Also, the presence of the heat transition and ice line 
traps may be the origin of massive planets with (very) large 
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orbital radii around massive stars which are detected by the 
direct imaging method. 

(iii) We have demonstrated that these planet traps can 
interact with each other, as they follow the time-dependent, 
viscous evolution of the disc (see Fig. (6]). More specifically, 
we have found that their characteristic radii evolve with dif- 
ferent dependence on M as the disc accretion rate M drops. 
This result shows that different planet traps move inwards 
differently, as discs accrete onto the central stars and star 
formation is terminated. (Equivalently, this evolution results 
in decreasing the surface density of gas (So) and accretion 
rates (M).) Our analyses are important for comprehensively 
understanding the formation of planetary systems around 
stars with various masses (see Fig. [7]). Since trapped planets 
can grow quickly, the interaction of active traps can initiate 
planet-planet interaction, which is controlled by the most 
massive, trapped planet. Thus, it is important to investigate 
how trapped planets gain their masses in viscously evolving 
discs. We will address this issue in a forthcoming paper. 

(iv) We have applied our results to planet traps in discs 
around massive (Herbig Ae/Be), intermediate (CTTS), and 
low (M) mass stars. In all cases, the heat transition trap is 
the outermost trap in the disc while the dead zone or the 
ice line trap is the innermost. At low accretion rates (M ~ 
10~^ Mq /year) in CTTS system as an example, the traps 
are located at r^dge ~ 1.2 au, ru — 1.5 au, and rht ~ 2.4 au, 
respectively (see Table [9|. 

(v) We have shown that the position of planet traps 
strongly depends on stellar masses and disc accretion rates. 
This indicates that host stars dictate a preferred scale of 
planetary systems formed around them (see Table [9|. For 
Iferbig Ae/Be stars, the size of planetary systems may ex- 
pand from the order of a few au to ~ 10'^ au. For CTTSs, it 
shrinks to a few ten au. For M stars, any planetary systems 
formed around them may be well confined within the order 
of au. 

(vi) We have applied our analyses to the Nice model for 
Solar system evolution. By choosing the appropriate param- 
eters {sA,t, Sao), we have shown that the ice line trap can 
capture the core of Jupiter while the dead zone trap Saturn. 
The locations of these trapped cores match those of the Nice 
model at M = 5.5 x W~^Mq /year which is achieved towards 
the end of disc evolution. Furthermore, our analyses suggest 
why Jupiter is more massive than Saturn, which is one of 
the most crucial assumptions adopted in the improved Nice 
model. 

In a subsequent paper, we will investigate the growth 
of planets which are captured at the planet traps during the 
time-dependent, viscous evolution of the disc. Inclusion of 
planetary growth will allow us to address when the planet- 
planet interaction dominates planetary migration. 
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APPENDIX A: EXTENSION OF 
PAARDEKOOPER'S TORQUE FORMULA 

Here, we discuss the extension of Paardekooper's torque 
formula that is valid exclusively in pure power-law discs 
(see equation p.2p ). Since entropy gradients that scale an 
entropy-related corotation torque can be straightforwardly 
obtained from the disc profiles, we focus on the extension 
of Lindbald and vortensity-related corotation torques. For 
disc models, we adopt the same one in § 13.11 that is, gen- 
eral profiles for the surface density (see equation (|3.ip ) and 
power-law profiles for the disc temperature (T oc r'). We 
only examine the case of s = — 1 for Ei„t (see equation p.lf) ). 



because the results for the MMSN case are qualitatively sim- 
ilar. We also discuss the effects of the disc temperature that 
has more general profiles. 

For the Lindblad torque, we indeed derive an analytical 
relation that controls the direction of migration for the above 
disc model in § 14.31 (see equation (|4.20[) '). Exploiting the 
analysis, we estimate the critical value of Fg^crit above which 
the Lindblad torque dictates outward migration: 

_ tanh(l/e) + /ip(g-t/2 + 7/4) 

tanh(l/c)-/ip(s- t/2 + 7/4)' ^ ' 

where hp is the disc aspect ratio (H/r) at the position of a 
planet (see l4.3.2l for a full derivation). Adopting hp = 0.1, 
we find that Fg^crit = 1-5 for c = 1 while Fg^crit = 2.7 for 
c = 3. As discussed in § 14.3.21 the direction of migration is 
determined only by the Lindblad torque if the density mod- 
ification is larger than Fg^crit- This is because the Lindblad 
resonances distribute much further way from a planet than 
the horseshoe region. Based on the results of MG04 (see their 
fig. 1), the density distortion that is produced by the opacity 
and heat transitions is likely less than Fg^crit- In the follow- 
ing discussion, we adopt these critical values of Fg^crit- Thus, 
Lindblad torques surely become weaker around the opacity 
and heat transitions, but never give the dominant contri- 
bution to the direction of migration. Therefore, we ad opt 
the original term derived in [Paardekooper et al.l (|201G| ) for 
consistency in the usage of the torque formulation. 

For the vortensity-related corotation torque, we add a 
correction factor in the original formulation (see equa- 
tions (|3.5|l and (|3.3p ). This is because gradients of vortensity 
that are the core of the vortensity-related corotation torque, 
are very sensitive to the disc structure. Fig. lAll shows the 
surface density structure, the resultant behaviour of , and 
the critical value of t below which planets migrate outwards 
(see equation (|3.5p ) on the top, middle, and bottom panel, 
respectively. Since the vortensity gradients depend on the 
disc temperature through fl (see equations (|3.4|) and (|4.3p ). 
self-consistent treatments are required to derive t. For sim- 
plicity, however, we fix t = —1.2 for determining 0^ in this 
figure. This profile is established for viscously heated, opti- 
cally thick discs (see § I3.2.2|l . We also discuss the effects of 
the disc temperature with more general profiles below. 

Fig. lAll shows very complicated behaviour for (f)v ■ It is 
obvious that the choice of c = 1 (the solid line) is more 
appropriate for the opacity and heat transitions based on 
the results of MG04. However, we also consider c = 3 (the 
dashed line) for a parameter study. By comparing the mid- 
dle and bottom panels, one can observe immediately that 
the crucial effects of cpv on t arise from Gaussian-like func- 
tions that distribute around r trans- It is important that the 
extent of the functions is well characterised by the transi- 
tion width (oj = cH) (see the vertical lines on the bottom 
panel). For the case of c = 3, the function expands beyond 
the transition width. However, this is a result of our choice 
of the temperature profile. Adopting more general profiles 
similar to equation p.ip for the disc temperature, the ex- 
tent of this function strongly diminishes (see Fig. lA2p . Thus, 
it is reasonable to conclude that the effects of are very 
local and well confined with the transition width. Another 
important point is that the required value of t is strongly re- 
duced by <;/>„ . Such small values of t are never accomplished in 
protoplanetary discs. This means that the vortensity-related 
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corotation torque becomes significantly large and results in 
inward migration. Therefore, inclusion of 4>v reduces the pos- 
sible outward migration region which is a consequence of 
viscous heating and entropy-related corotation torques. 

In summary, proper treatments of the vortensity-related 
corotation torque are important for discs with general pro- 
files. However, the deviation from the power-law treatments 
is well confined in the vicinity of the position of disc inhomo- 
geneities, and therefore power-law approximation provides 
reliable results in our level of extension. 



APPENDIX B: TAYLOR EXPANSION OF THE 
TIDAL TORQUE 

Here, we derive the Taylor expansion of quantities related 
with the tidal torque in terms of 1/m <^ 1 since the torque 
takes the maximum value at the wavenumber m ~ 10 (IWardI 
1 19971 ). 
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Taking the limit ^ — >■ (since the effect of gas pressure 
becomes negligible in the Keplerian discs). 
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APPENDIX C: ICE LINES OF OTHER 
MOLECULES 

Here, we investigate the minimum abundance of a molecule 
that is required for its ice line to be a barrier, by using the 
analyses done in ij |4.4l Assuming the dust density to increase 
by the amount of ice which is formed at the ice line radius of 
a specific molecule, r, the required dust density enhancement 
is (using equations (|4.35|l and (|4.27|) ). for r/r^dge <1, 



fdi > 



= f{aA,aD,SA,t,h,r,fd2,c), 



(CI) 




2.0 



Figure Al. The structure of the vortensity correction factor (j>y 
for discs with general profiles. On the top panel, the surface den- 
sity structures are plotted (see equation l|3.1| lV On the middle, the 
resultant behaviours of (py are shown. On the bottom, the critical 
values of t affected by are plotted. We stress that the choice 
of c = 1 is more appropriate for opacity and heat transitions. For 
the purpose of a parameter study, we also consider c = 3. The 
results of pure power-law discs are denoted by the horizontal lines 
on the middle and bottom panels. The crucial effects of <j) y arise 
from Gaussian-like functions that are well confined with the tran- 
sition width {uj = cH) (also see Fig. |A2j. Inclusion of (pv signifi- 
cantly lowers the required value of t around r = rtrans- This indi- 
cates that the vortensity-related corotation torque becomes con- 
siderably larger there and makes planets migrate inwards. Thus, 
proper treatments of the vortensity-related corotation torque de- 
viates the results of the pure power-law discs, but the deviations 
occur only in the local region centered at r = rtrans- 
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Figure A2. The critical value of t for the disc temperature with 
general profiles. We adopt an expression similar to equation Il3.1|l 
for the disc temperature. Also, we consider the case of c = 1. By 
comparing Fig. lAll the extent wherein 4iv plays a crucial role is 
strongly diminished. This supports the usage of the pure power- 
law approximation. 
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For rjr^dge >1, the dust density enhancement needed is 
given as 

-3t/2 + 1/4 



9d > 



(C5) 



3t/2 - 1/4 + Vhi2/c/i 

where equation (|4.4H) is used. In order to gain a unified 
condition for the entire of discs, we find the relation between 
fdi and Qd which is given as 
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As a result, equation (|C5|I becomes 
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Fig. ICll shows the threshold value of fdi as a 
function of r/r^dge- We define a fiducial model as 
(sA,t, /ledge, /d2,c)=(3,-l/2,0. 1,0,1). The required value of 
fdi are above the lines. We performed parameter studies on 
hedge, SA, and t and the results are shown in the top, middle, 
and bottom panels, respectively. 

For r/r-edge <1, the dominant term in equation HC1|I is 
[redge/ry^^^^^^"^ in Kr. We see that a large value of fdi is 
required for outward migration, if the value of SA + t + 3/2 is 
positive. Fig. ICll shows that small separations between the 



ice line and dead zone barriers are again preferred. If SA + i + 
3/2 <0, a tiny amount of ice is enough to provide a barrier 
(see the dashed line on the middle panel). This indicates 
that ice lines of relatively abundant molecules can work as 
barriers. For r/redge >1, the dominant term in equation 
(fC7)l is {r/rf.dge)"^^'^*^^^ in Pi. Consequently, a significantly 
large value of fdi is needed for outward migration, if the 
value oi SA + t -\- 3/2 is positive. In fact, it may not be 
possible for ice lines of any molecule to be a barrier in this 
case. On the other hand, if sa + t + 3/2 is negative, the 
amount of ice does not matter, (see the dashed line on the 
middle panel). This implies that ice lines of any species can 
act as barriers. Thus, the required fdi strongly depends on 
the sum of sa and t. 

What value of sa is most appropriate for realistic discs? 
For r/redge <1, Sa should be positive. Otherwise, dead zones 
cannot exist. This indicates that larger values of fdi are re- 
quired as r/r^dge decreases, and hence only specific, abun- 
dant species in discs can work as a barrier. For r/redge >1, 
Sa should be negative, since Ea ~ E. Thus, it is likely that 
the abundance of molecules may not be important for pro- 
viding a barrier. As a result, several ice line barriers which 
are excited by different species such as water-ice and CO-ice 
may be present in a single disc. However, the value of sa 
strongly depends on how icy dust absorbs free electrons and 
how these distributions evolve with time. 

This paper has been typeset from a Tj^X/ I^Tj^X file prepared 
by the author. 
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Figure CI. The enhancement f^i of the dust density required 
to produce ice Hne barriers. For r/r^^g^ <1, equation lICll l is 
adopted while equation ||C7| | is used for r/r^^g^ >1. We define 
a fiducial model as (syi, t, /ij^jge, /j2i c)=(3,-l/2,0.1,0,l). We vary 
hedge, ^A, and t On the top, middle, and bottom panel, respec- 
tively. A large value of /^ji is required for ice lines to be a barrier, 
if Syl + 1 + 3/2 is positive. On the other hand, the presence of ice 
line barriers becomes insensitive to f^i if + * + 3/2 is nega- 
tive (see the dashed line on the middle panel). It implies that the 
minimum value of f^i stronelv depends on the detail structure of 



